1 /* @source primersearch application
2 **
3 ** Searches a set of primer pairs against a set of DNA sequences in both
4 ** forward and reverse sense.
5 ** Modification of fuzznuc.
6 ** @author Copyright (C) Val Curwen (vac@sanger.ac.uk)
7 ** @@
8 **
9 ** This program is free software; you can redistribute it and/or
10 ** modify it under the terms of the GNU General Public License
11 ** as published by the Free Software Foundation; either version 2
12 ** of the License, or (at your option) any later version.
13 **
14 ** This program is distributed in the hope that it will be useful,
15 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
16 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 ** GNU General Public License for more details.
18 **
19 ** You should have received a copy of the GNU General Public License
20 ** along with this program; if not, write to the Free Software
21 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
22 ******************************************************************************/
23
24
25 #include "emboss.h"
26 #include <stdlib.h>
27 #include <string.h>
28 #include <stdio.h>
29
30 #define AJA2 50
31
32
33
34
35 /* @datastatic PGuts **********************************************************
36 **
37 ** the internals of a primer; each Primer has two of these,
38 ** one forward and one reverse
39 **
40 ** @alias primerguts
41 **
42 ** @attr patstr [AjPStr] Undocumented
43 ** @attr origpat [AjPStr] Undocumented
44 ** @attr type [ajuint] Undocumented
45 ** @attr len [ajuint] Undocumented
46 ** @attr real_len [ajuint] Undocumented
47 ** @attr amino [AjBool] Undocumented
48 ** @attr carboxyl [AjBool] Undocumented
49 ** @attr mm [ajuint] Undocumented
50 ** @attr buf [ajint*] Undocumented
51 ** @attr sotable [ajuint*] Undocumented
52 ** @attr off [EmbOPatBYPNode[AJALPHA]] Undocumented
53 ** @attr re [AjPStr] Undocumented
54 ** @attr skipm [ajuint**] Undocumented
55 ** @attr tidy [const void*] Undocumented
56 ** @attr solimit [ajuint] Undocumented
57 ** @attr Padding [char[4]] Padding to alignment boundary
58 ******************************************************************************/
59
60 typedef struct primerguts
61 {
62 AjPStr patstr;
63 AjPStr origpat;
64 ajuint type;
65 ajuint len;
66 ajuint real_len;
67 AjBool amino;
68 AjBool carboxyl;
69
70 ajuint mm;
71
72 ajint* buf;
73 ajuint* sotable;
74
75 EmbOPatBYPNode off[AJALPHA];
76 AjPStr re;
77 ajuint** skipm;
78 const void* tidy;
79 ajuint solimit;
80 char Padding[4];
81 } *PGuts;
82
83
84
85
86 /* @datastatic PHit ***********************************************************
87 **
88 ** holds details of a hit against a sequence ie this primer will amplify
89 **
90 ** @alias primerhit
91 **
92 ** @attr seqname [AjPStr] Undocumented
93 ** @attr desc [AjPStr] Undocumented
94 ** @attr acc [AjPStr] Undocumented
95 ** @attr forward [AjPStr] pattern that hits forward strand
96 ** @attr reverse [AjPStr] pattern that hits reverse strand
97 ** @attr forward_pos [ajuint] Undocumented
98 ** @attr reverse_pos [ajuint] Undocumented
99 ** @attr amplen [ajuint] Undocumented
100 ** @attr forward_mismatch [ajuint] Undocumented
101 ** @attr reverse_mismatch [ajuint] Undocumented
102 ** @attr Padding [char[4]] Padding to alignment boundary
103 ******************************************************************************/
104
105 typedef struct primerhit
106 {
107 AjPStr seqname;
108 AjPStr desc;
109 AjPStr acc;
110 AjPStr forward;
111 AjPStr reverse;
112 ajuint forward_pos;
113 ajuint reverse_pos;
114 ajuint amplen;
115 ajuint forward_mismatch;
116 ajuint reverse_mismatch;
117 char Padding[4];
118 } *PHit;
119
120
121
122
123 /* @datastatic Primer *********************************************************
124 **
125 ** primer pairs will be read into a list of these structs
126 **
127 ** @alias primers
128 **
129 ** @attr Name [AjPStr] Undocumented
130 ** @attr forward [PGuts] Undocumented
131 ** @attr reverse [PGuts] Undocumented
132 ** @attr hitlist [AjPList] Undocumented
133 ******************************************************************************/
134
135 typedef struct primers
136 {
137 AjPStr Name;
138 PGuts forward;
139 PGuts reverse;
140 AjPList hitlist;
141 } *Primer;
142
143
144
145
146 /* "constructors" */
147 static void primersearch_initialise_pguts(PGuts* primer);
148
149 /* "destructors" */
150 static void primersearch_free_pguts(PGuts* primer);
151 static void primersearch_free_primer(void** x, void* cl);
152 static void primersearch_clean_hitlist(AjPList* hlist);
153
154 /* utilities */
155 static void primersearch_read_primers(AjPList* primerList, AjPFile primerFile,
156 ajint mmp);
157 static AjBool primersearch_classify_and_compile(Primer* primdata);
158 static void primersearch_primer_search(const AjPList primerList,
159 const AjPSeq seq);
160 static void primersearch_scan_seq(const Primer primdata,
161 const AjPSeq seq, AjBool reverse);
162 static void primersearch_store_hits(const Primer primdata, AjPList fhits_list,
163 AjPList rhits_list, const AjPSeq seq,
164 AjBool reverse);
165 static void primersearch_print_hits(const AjPList primerList, AjPFile outf);
166
167
168
169
170 /* @prog primersearch *********************************************************
171 **
172 ** Searches DNA sequences for matches with primer pairs
173 **
174 ******************************************************************************/
175
main(int argc,char ** argv)176 int main(int argc, char **argv)
177 {
178 AjPSeqall seqall;
179 AjPSeq seq = NULL;
180 AjPFile primerFile; /* read the primer pairs from a file */
181 AjPFile outf;
182 AjPList primerList;
183
184 ajint mmp = 0;
185
186 embInit("primersearch", argc, argv);
187
188 seqall = ajAcdGetSeqall("seqall");
189 outf = ajAcdGetOutfile("outfile");
190 primerFile = ajAcdGetInfile("infile");
191 mmp = ajAcdGetInt("mismatchpercent");
192
193 /* build list of forward/reverse primer pairs as read from primerfile */
194 primerList = ajListNew();
195
196 /* read in primers from primerfile, classify and compile them */
197 primersearch_read_primers(&primerList,primerFile, mmp);
198
199 /* check there are primers to be searched */
200 if(!ajListGetLength(primerList))
201 {
202 ajErr("No suitable primers found - exiting");
203 embExitBad();
204 return 0;
205
206 }
207
208 /* query sequences one by one */
209 while(ajSeqallNext(seqall,&seq))
210 primersearch_primer_search(primerList, seq);
211
212 /* output the results */
213 primersearch_print_hits(primerList, outf);
214
215 /* delete all nodes of list, then the list itself */
216 ajListMap(primerList, &primersearch_free_primer, NULL);
217 ajListFree(&primerList);
218 ajListFree(&primerList);
219
220 ajFileClose(&outf);
221
222 ajSeqallDel(&seqall);
223 ajSeqDel(&seq);
224
225 ajFileClose(&primerFile);
226
227 embExit();
228
229 return 0;
230 }
231
232
233
234
235 /* "constructors" */
236
237
238
239
240 /* @funcstatic primersearch_initialise_pguts **********************************
241 **
242 ** Initialise primer guts
243 **
244 ** @param [w] primer [PGuts*] primer guts
245 ** @@
246 ******************************************************************************/
247
primersearch_initialise_pguts(PGuts * primer)248 static void primersearch_initialise_pguts(PGuts* primer)
249 {
250
251 AJNEW(*primer);
252 (*primer)->patstr = NULL;
253 (*primer)->origpat = ajStrNew();
254 (*primer)->type = 0;
255 (*primer)->len = 0;
256 (*primer)->real_len = 0;
257 (*primer)->re = NULL;
258 (*primer)->amino = 0;
259 (*primer)->carboxyl = 0;
260 (*primer)->tidy = NULL;
261
262 (*primer)->mm = 0;
263 (*primer)->buf = NULL;
264 (*primer)->sotable = NULL;
265 (*primer)->solimit = 0;
266 (*primer)->re = NULL;
267 (*primer)->skipm = NULL;
268
269 return;
270 }
271
272
273
274
275 /* "destructors" */
276
277
278
279
280 /* @funcstatic primersearch_free_pguts ****************************************
281 **
282 ** Frees up all the internal members of a PGuts struct
283 **
284 ** @param [d] primer [PGuts*] Undocumented
285 ** @return [void]
286 ** @@
287 ******************************************************************************/
288
primersearch_free_pguts(PGuts * primer)289 static void primersearch_free_pguts(PGuts* primer)
290 {
291 ajuint i = 0;
292
293 ajStrDel(&(*primer)->patstr);
294 ajStrDel(&(*primer)->origpat);
295 ajStrDel(&(*primer)->re);
296
297
298 if(((*primer)->type==1 || (*primer)->type==2) && ((*primer)->buf))
299 free((*primer)->buf);
300
301 if(((*primer)->type==3 || (*primer)->type==4) && ((*primer)->sotable))
302 free((*primer)->sotable);
303
304 if((*primer)->type==6)
305 for(i=0;i<(*primer)->real_len;++i) AJFREE((*primer)->skipm[i]);
306 AJFREE(*primer);
307
308 return;
309 }
310
311
312
313
314 /* @funcstatic primersearch_free_primer ***************************************
315 **
316 ** frees up the internal members of a Primer
317 **
318 ** @param [r] x [void**] Undocumented
319 ** @param [r] cl [void*] Undocumented
320 ** @@
321 ******************************************************************************/
322
primersearch_free_primer(void ** x,void * cl)323 static void primersearch_free_primer(void **x, void *cl)
324 {
325 Primer* p;
326 Primer primdata;
327 AjIList lIter;
328
329 (void) cl; /* make it used */
330
331 p = (Primer*) x;
332 primdata = *p;
333
334 primersearch_free_pguts(&primdata->forward);
335 primersearch_free_pguts(&primdata->reverse);
336 ajStrDel(&primdata->Name);
337
338 /* clean up hitlist */
339 lIter = ajListIterNewread(primdata->hitlist);
340 while(!ajListIterDone(lIter))
341 {
342 PHit phit = ajListIterGet(lIter);
343 ajStrDel(&phit->forward);
344 ajStrDel(&phit->reverse);
345 ajStrDel(&phit->seqname);
346 ajStrDel(&phit->acc);
347 ajStrDel(&phit->desc);
348
349 AJFREE(phit);
350 }
351
352 ajListFree(&primdata->hitlist);
353 ajListFree(&primdata->hitlist);
354 ajListIterDel(&lIter);
355
356 AJFREE(primdata);
357
358 return;
359 }
360
361
362
363
364 /* @funcstatic primersearch_clean_hitlist *************************************
365 **
366 ** Clean the hitlist
367 **
368 ** @param [d] hlist [AjPList*] Undocumented
369 ** @@
370 ******************************************************************************/
371
primersearch_clean_hitlist(AjPList * hlist)372 static void primersearch_clean_hitlist(AjPList* hlist)
373 {
374 AjIList lIter;
375
376 lIter = ajListIterNewread(*hlist);
377 while(!ajListIterDone(lIter))
378 {
379 EmbPMatMatch fm = ajListIterGet(lIter);
380 embMatMatchDel(&fm);
381 }
382 ajListFree(hlist);
383 ajListFree(hlist);
384 ajListIterDel(&lIter);
385
386 return;
387 }
388
389
390
391
392 /* utilities */
393
394
395
396
397 /* @funcstatic primersearch_read_primers **************************************
398 **
399 ** read primers in from primerfile, classify and compile the patterns
400 **
401 ** @param [w] primerList [AjPList*] primer list
402 ** @param [u] primerFile [AjPFile] Undocumented
403 ** @param [r] mmp [ajint] Undocumented
404 ** @@
405 ******************************************************************************/
406
primersearch_read_primers(AjPList * primerList,AjPFile primerFile,ajint mmp)407 static void primersearch_read_primers(AjPList *primerList, AjPFile primerFile,
408 ajint mmp)
409 {
410 AjPStr rdline = NULL;
411 AjPStrTok handle = NULL;
412
413 ajint nprimers = 0;
414 Primer primdata = NULL;
415
416
417 while (ajReadlineTrim(primerFile, &rdline))
418 {
419 primdata = NULL;
420 if (ajStrGetCharFirst(rdline) == '#')
421 continue;
422 if (ajStrSuffixC(rdline, ".."))
423 continue;
424
425 AJNEW(primdata);
426 primdata->Name = NULL;
427
428 primersearch_initialise_pguts(&primdata->forward);
429 primersearch_initialise_pguts(&primdata->reverse);
430
431 primdata->hitlist = ajListNew();
432
433 handle = ajStrTokenNewC(rdline, " \t");
434 ajStrTokenNextParse(handle, &primdata->Name);
435
436 ajStrTokenNextParse(handle, &primdata->forward->patstr);
437 ajStrFmtUpper(&primdata->forward->patstr);
438 ajStrTokenNextParse(handle, &primdata->reverse->patstr);
439 ajStrFmtUpper(&primdata->reverse->patstr);
440 ajStrTokenDel(&handle);
441
442 /* copy patterns for Henry Spencer code */
443 ajStrAssignC(&primdata->forward->origpat,
444 ajStrGetPtr(primdata->forward->patstr));
445 ajStrAssignC(&primdata->reverse->origpat,
446 ajStrGetPtr(primdata->reverse->patstr));
447
448 /* set the mismatch level */
449 primdata->forward->mm = (ajint)
450 (ajStrGetLen(primdata->forward->patstr)*mmp)/100;
451 primdata->reverse->mm = (ajint)
452 (ajStrGetLen(primdata->reverse->patstr)*mmp)/100;
453
454 if(primersearch_classify_and_compile(&primdata))
455 {
456 ajListPushAppend(*primerList, primdata);
457 nprimers++;
458 }
459 else /* there was something funny about the primer sequences */
460 {
461 ajUser("Cannot use %s\n", ajStrGetPtr(primdata->Name));
462 primersearch_free_pguts(&primdata->forward);
463 primersearch_free_pguts(&primdata->reverse);
464 ajStrDel(&primdata->Name);
465 ajListFree(&primdata->hitlist);
466 ajListFree(&primdata->hitlist);
467 AJFREE(primdata);
468 }
469 }
470
471 ajStrDel(&rdline);
472
473 return;
474 }
475
476
477
478
479 /* @funcstatic primersearch_classify_and_compile ******************************
480 **
481 ** determines pattern type and compiles it
482 **
483 ** @param [w] primdata [Primer*] primer data
484 ** @return [AjBool] true if useable primer
485 ** @@
486 ******************************************************************************/
487
primersearch_classify_and_compile(Primer * primdata)488 static AjBool primersearch_classify_and_compile(Primer* primdata)
489 {
490
491 /* forward primer */
492 if(!((*primdata)->forward->type =
493 embPatGetType(((*primdata)->forward->origpat),
494 &((*primdata)->forward->patstr),
495 (*primdata)->forward->mm,0,
496 &((*primdata)->forward->real_len),
497 &((*primdata)->forward->amino),
498 &((*primdata)->forward->carboxyl))))
499 ajFatal("Illegal pattern");
500
501 /* reverse primer */
502 if(!((*primdata)->reverse->type =
503 embPatGetType(((*primdata)->reverse->origpat),
504 &((*primdata)->reverse->patstr),
505 (*primdata)->reverse->mm,0,
506 &((*primdata)->reverse->real_len),
507 &((*primdata)->reverse->amino),
508 &((*primdata)->reverse->carboxyl))))
509 ajFatal("Illegal pattern");
510
511 embPatCompile((*primdata)->forward->type,
512 (*primdata)->forward->patstr,
513 &((*primdata)->forward->len),
514 &((*primdata)->forward->buf),
515 (*primdata)->forward->off,
516 &((*primdata)->forward->sotable),
517 &((*primdata)->forward->solimit),
518 &((*primdata)->forward->real_len),
519 &((*primdata)->forward->re),
520 &((*primdata)->forward->skipm),
521 (*primdata)->forward->mm );
522
523 embPatCompile((*primdata)->reverse->type,
524 (*primdata)->reverse->patstr,
525 &((*primdata)->reverse->len),
526 &((*primdata)->reverse->buf),
527 (*primdata)->reverse->off,
528 &((*primdata)->reverse->sotable),
529 &((*primdata)->reverse->solimit),
530 &((*primdata)->reverse->real_len),
531 &((*primdata)->reverse->re),
532 &((*primdata)->reverse->skipm),
533 (*primdata)->reverse->mm );
534
535 return AJTRUE; /* this is a useable primer */
536 }
537
538
539
540
541 /* @funcstatic primersearch_primer_search *************************************
542 **
543 ** tests the primers in primdata against seq and writes results to outfile
544 **
545 ** @param [r] primerList [const AjPList] primer list
546 ** @param [r] seq [const AjPSeq] sequence
547 ** @@
548 ******************************************************************************/
549
primersearch_primer_search(const AjPList primerList,const AjPSeq seq)550 static void primersearch_primer_search(const AjPList primerList,
551 const AjPSeq seq)
552 {
553 AjIList listIter;
554
555 /* test each list node against this sequence */
556 listIter = ajListIterNewread(primerList);
557 while(!ajListIterDone(listIter))
558 {
559 Primer curr_primer = ajListIterGet(listIter);
560
561 primersearch_scan_seq(curr_primer, seq, AJFALSE);
562 primersearch_scan_seq(curr_primer, seq, AJTRUE);
563 }
564
565 ajListIterDel(&listIter);
566
567 return;
568 }
569
570
571
572
573 /* @funcstatic primersearch_scan_seq ******************************************
574 **
575 ** scans the primer pairs against the sequence in either forward
576 ** sense or reverse complemented
577 ** works out amplimer length if the two primers both hit
578 **
579 ** @param [r] primdata [const Primer] primer data
580 ** @param [r] seq [const AjPSeq] sequence
581 ** @param [r] reverse [AjBool] do reverse
582 ** @@
583 ******************************************************************************/
584
primersearch_scan_seq(const Primer primdata,const AjPSeq seq,AjBool reverse)585 static void primersearch_scan_seq(const Primer primdata,
586 const AjPSeq seq, AjBool reverse)
587 {
588 AjPStr seqstr = NULL;
589 AjPStr revstr = NULL;
590 AjPStr seqname = NULL;
591 ajuint fhits = 0;
592 ajuint rhits = 0;
593 AjPList fhits_list = NULL;
594 AjPList rhits_list = NULL;
595
596 /* initialise variables for search */
597 ajStrAssignC(&seqname,ajSeqGetNameC(seq));
598 ajStrAssignS(&seqstr, ajSeqGetSeqS(seq));
599 ajStrAssignS(&revstr, ajSeqGetSeqS(seq));
600 ajStrFmtUpper(&seqstr);
601 ajStrFmtUpper(&revstr);
602 ajSeqstrReverse(&revstr);
603 fhits_list = ajListNew();
604 rhits_list = ajListNew();
605
606 if(!reverse)
607 {
608 /* test OligoA against forward sequence, and OligoB against reverse */
609 embPatFuzzSearch(primdata->forward->type,
610 ajSeqGetBegin(seq),
611 primdata->forward->patstr,
612 seqname,
613 seqstr,
614 fhits_list,
615 primdata->forward->len,
616 primdata->forward->mm,
617 primdata->forward->amino,
618 primdata->forward->carboxyl,
619 primdata->forward->buf,
620 primdata->forward->off,
621 primdata->forward->sotable,
622 primdata->forward->solimit,
623 primdata->forward->re,
624 primdata->forward->skipm,
625 &fhits,
626 primdata->forward->real_len,
627 &(primdata->forward->tidy));
628
629 if(fhits)
630 embPatFuzzSearch(primdata->reverse->type,
631 ajSeqGetBegin(seq),
632 primdata->reverse->patstr,
633 seqname,
634 revstr,
635 rhits_list,
636 primdata->reverse->len,
637 primdata->reverse->mm,
638 primdata->reverse->amino,
639 primdata->reverse->carboxyl,
640 primdata->reverse->buf,
641 primdata->reverse->off,
642 primdata->reverse->sotable,
643 primdata->reverse->solimit,
644 primdata->reverse->re,
645 primdata->reverse->skipm,
646 &rhits,
647 primdata->reverse->real_len,
648 &(primdata->reverse->tidy));
649 }
650 else
651 {
652 /*test OligoB against forward sequence, and OligoA against reverse */
653 embPatFuzzSearch(primdata->reverse->type,
654 ajSeqGetBegin(seq),
655 primdata->reverse->patstr,
656 seqname,
657 seqstr,
658 fhits_list,
659 primdata->reverse->len,
660 primdata->reverse->mm,
661 primdata->reverse->amino,
662 primdata->reverse->carboxyl,
663 primdata->reverse->buf,
664 primdata->reverse->off,
665 primdata->reverse->sotable,
666 primdata->reverse->solimit,
667 primdata->reverse->re,
668 primdata->reverse->skipm,
669 &fhits,
670 primdata->reverse->real_len,
671 &(primdata->reverse->tidy));
672
673 if(fhits)
674 embPatFuzzSearch(primdata->forward->type,
675 ajSeqGetBegin(seq),
676 primdata->forward->patstr,
677 seqname,
678 revstr,
679 rhits_list,
680 primdata->forward->len,
681 primdata->forward->mm,
682 primdata->forward->amino,
683 primdata->forward->carboxyl,
684 primdata->forward->buf,
685 primdata->forward->off,
686 primdata->forward->sotable,
687 primdata->forward->solimit,
688 primdata->forward->re,
689 primdata->forward->skipm,
690 &rhits,
691 primdata->forward->real_len,
692 &(primdata->forward->tidy));
693 }
694
695 if(fhits && rhits)
696 /* get amplimer length(s) and write out the hit */
697 primersearch_store_hits(primdata, fhits_list, rhits_list,
698 seq, reverse);
699
700 /* tidy up */
701 primersearch_clean_hitlist(&fhits_list);
702 primersearch_clean_hitlist(&rhits_list);
703
704 ajStrDel(&seqstr);
705 ajStrDel(&revstr);
706 ajStrDel(&seqname);
707
708 return;
709 }
710
711
712
713
714 /* @funcstatic primersearch_store_hits ****************************************
715 **
716 ** Store primer hits
717 **
718 ** @param [r] primdata [const Primer] primer data
719 ** @param [w] fhits [AjPList] forward hits
720 ** @param [w] rhits [AjPList] reverse hits
721 ** @param [r] seq [const AjPSeq] sequence
722 ** @param [r] reverse [AjBool] do reverse
723 ** @@
724 ******************************************************************************/
725
primersearch_store_hits(const Primer primdata,AjPList fhits,AjPList rhits,const AjPSeq seq,AjBool reverse)726 static void primersearch_store_hits(const Primer primdata,
727 AjPList fhits, AjPList rhits,
728 const AjPSeq seq, AjBool reverse)
729 {
730 ajint amplen = 0;
731 AjIList fi;
732 AjIList ri;
733
734 PHit primerhit = NULL;
735
736 fi = ajListIterNewread(fhits);
737 while(!ajListIterDone(fi))
738 {
739 EmbPMatMatch fm = NULL;
740 EmbPMatMatch rm = NULL;
741 amplen = 0;
742
743 fm = ajListIterGet(fi);
744 ri = ajListIterNewread(rhits);
745 while(!ajListIterDone(ri))
746 {
747 ajint seqlen = ajSeqGetLen(seq);
748 ajint s = (fm->start);
749 ajint e;
750
751 rm = ajListIterGet(ri);
752 e = (rm->start-1);
753 amplen = seqlen-(s-1)-e;
754
755 if (amplen > 0) /* no point making a hit if -ve length! */
756 {
757 primerhit = NULL;
758 AJNEW(primerhit);
759 primerhit->desc=NULL; /* must be NULL for ajStrAss */
760 primerhit->seqname=NULL; /* must be NULL for ajStrAss */
761 primerhit->acc=NULL;
762 primerhit->forward=NULL;
763 primerhit->reverse=NULL;
764 ajStrAssignC(&primerhit->seqname,ajSeqGetNameC(seq));
765 ajStrAssignS(&primerhit->desc, ajSeqGetDescS(seq));
766 ajStrAssignS(&primerhit->acc, ajSeqGetAccS(seq));
767 primerhit->forward_pos = fm->start;
768 primerhit->reverse_pos = rm->start;
769 primerhit->forward_mismatch = fm->mm;
770 primerhit->reverse_mismatch = rm->mm;
771 primerhit->amplen = amplen;
772 if(!reverse)
773 {
774 ajStrAssignS(&primerhit->forward,
775 primdata->forward->patstr);
776 ajStrAssignS(&primerhit->reverse,
777 primdata->reverse->patstr);
778 }
779 else
780 {
781 ajStrAssignS(&primerhit->forward,
782 primdata->reverse->patstr);
783 ajStrAssignS(&primerhit->reverse,
784 primdata->forward->patstr);
785 }
786 ajListPushAppend(primdata->hitlist, primerhit);
787
788
789 }
790 }
791 /*
792 ** clean up rListIter here as it will be new'ed again next
793 ** time through
794 */
795 ajListIterDel(&ri);
796 }
797
798 ajListIterDel(&fi);
799 return;
800 }
801
802
803
804
805 /* @funcstatic primersearch_print_hits ****************************************
806 **
807 ** Print primer hits
808 **
809 ** @param [r] primerList [const AjPList] primer hits
810 ** @param [w] outf [AjPFile] outfile
811 ** @@
812 ******************************************************************************/
813
primersearch_print_hits(const AjPList primerList,AjPFile outf)814 static void primersearch_print_hits(const AjPList primerList, AjPFile outf)
815 {
816 /* iterate through list of hits */
817 AjIList lIter;
818
819 ajint count = 1;
820 lIter = ajListIterNewread(primerList);
821 while(!ajListIterDone(lIter))
822 {
823 Primer primer = ajListIterGet(lIter);
824 AjIList hIter = ajListIterNewread(primer->hitlist);
825 count = 1;
826
827 ajFmtPrintF(outf, "\nPrimer name %s\n", ajStrGetPtr(primer->Name));
828
829 while(!ajListIterDone(hIter))
830 {
831 PHit hit = ajListIterGet(hIter);
832 ajFmtPrintF(outf, "Amplimer %d\n", count);
833 ajFmtPrintF(outf, "\tSequence: %s %s \n\t%s\n",
834 ajStrGetPtr(hit->seqname),
835 ajStrGetPtr(hit->acc), ajStrGetPtr(hit->desc));
836 ajFmtPrintF(outf, "\t%s hits forward strand at %d with %d "
837 "mismatches\n",
838 ajStrGetPtr(hit->forward), hit->forward_pos,
839 hit->forward_mismatch);
840 ajFmtPrintF(outf, "\t%s hits reverse strand at [%d] with %d "
841 "mismatches\n",
842 ajStrGetPtr(hit->reverse), (hit->reverse_pos),
843 (hit->reverse_mismatch));
844 ajFmtPrintF(outf, "\tAmplimer length: %d bp\n", hit->amplen);
845 count++;
846 }
847 ajListIterDel(&hIter);
848 }
849 ajListIterDel(&lIter);
850
851 return;
852 }
853