1 /* @source twofeat application
2 **
3 ** Finds neighbouring pairs of features in sequences
4 ** @author Copyright (C) Gary Williams (gwilliam@hgmp.mrc.ac.uk)
5 ** @@
6 **
7 ** This program is free software; you can redistribute it and/or
8 ** modify it under the terms of the GNU General Public License
9 ** as published by the Free Software Foundation; either version 2
10 ** of the License, or (at your option) any later version.
11 **
12 ** This program is distributed in the hope that it will be useful,
13 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
14 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 ** GNU General Public License for more details.
16 **
17 ** You should have received a copy of the GNU General Public License
18 ** along with this program; if not, write to the Free Software
19 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
20 ******************************************************************************/
21
22 #include "emboss.h"
23
24
25
26
27 /* @datastatic PHit ***********************************************************
28 **
29 ** twofeat internals
30 **
31 ** @alias SHit
32 ** @alias OHit
33 **
34 ** @attr gfA [AjPFeature] Undocumented
35 ** @attr gfB [AjPFeature] Undocumented
36 ** @attr Start [ajint] Undocumented
37 ** @attr End [ajint] Undocumented
38 ** @attr distance [ajint] Undocumented
39 ** @attr Padding [char[4]] Padding to alignment boundary
40 ******************************************************************************/
41
42 typedef struct SHit
43 {
44 AjPFeature gfA;
45 AjPFeature gfB;
46 ajint Start;
47 ajint End;
48 ajint distance;
49 char Padding[4];
50 } OHit;
51 #define PHit OHit*
52
53
54
55
56 static void twofeat_rippledown(const AjPFeattable tabA,
57 const AjPFeattable tabB,
58 ajint overlapi, ajint minrange,
59 ajint maxrange, ajint
60 rangetypei, ajint sensei,
61 AjPFeattable outtab,
62 AjBool twoout, const AjPStr typeout);
63 static AjBool twofeat_check_match(AjPFeature gfA, AjPFeature gfB,
64 PHit *detail, ajint overlapi,
65 ajint minrange, ajint maxrange, ajint
66 rangetypei, ajint sensei);
67 static void twofeat_report_hits(const AjPList hitlist, AjBool twoout,
68 const AjPStr typeout, AjPFeattable outtab);
69 static void twofeat_find_features(const AjPSeq seq, AjPFeattable tab,
70 ajint begin, ajint end, const AjPStr source,
71 const AjPStr type, ajint sense,
72 float minscore, float maxscore,
73 const AjPStr tag, const AjPStr value);
74 static AjBool twofeat_MatchFeature(const AjPFeature gf,
75 const AjPStr source, const AjPStr type,
76 ajint sense, float minscore,
77 float maxscore,
78 const AjPStr tag, const AjPStr value,
79 AjBool *tagsmatch);
80 static AjBool twofeat_MatchPatternTags(const AjPFeature feat,
81 const AjPStr tpattern,
82 const AjPStr vpattern);
83 static PHit twofeat_HitsNew(void);
84 static void twofeat_HitsDel(PHit *pthis);
85 static ajint twofeat_get_overlap_type(const AjPStr overlap);
86 static ajint twofeat_get_range_type(const AjPStr rangetype);
87 static ajint twofeat_get_sense_type(const AjPStr sense);
88 /* static ajint twofeat_get_order_type(const AjPStr order); */
89
90
91
92
93 /********* relation criterion types *****************/
94 enum OVERLAP_TYPE {OV_ANY, OV_OVERLAP, OV_NOTOVER, OV_NOTIN, OV_AIN, OV_BIN};
95 enum RANGE_TYPE {RA_NEAREST, RA_LEFT, RA_RIGHT, RA_FURTHEST};
96 enum SENSE_TYPE {SN_ANY, SN_SAME, SN_OPPOSITE};
97 enum ORDER_TYPE {OR_ANY, OR_AB, OR_BA};
98
99
100
101
102 /* @prog twofeat **************************************************************
103 **
104 ** Finds neighbouring pairs of features in sequences
105 **
106 ******************************************************************************/
107
main(int argc,char ** argv)108 int main(int argc, char **argv)
109 {
110 /* input */
111 AjPSeqall seqall;
112
113 /* feature a */
114 AjPStr asource;
115 AjPStr atype;
116 AjPStr asense;
117 ajint asensei;
118 ajint aminscore;
119 ajint amaxscore;
120 AjPStr atag;
121 AjPStr avalue;
122
123 /* feature b */
124 AjPStr bsource;
125 AjPStr btype;
126 AjPStr bsense;
127 ajint bsensei;
128 ajint bminscore;
129 ajint bmaxscore;
130 AjPStr btag;
131 AjPStr bvalue;
132
133 /* relation */
134 AjPStr overlap;
135 ajint minrange;
136 ajint maxrange;
137 AjPStr rangetype;
138 AjPStr sense;
139 AjPStr order;
140 ajint overlapi;
141 ajint rangetypei;
142 ajint sensei;
143
144 /* output */
145 AjBool twoout;
146 AjPStr typeout;
147 AjPReport report = NULL;
148
149
150 AjPSeq seq;
151 AjPFeattable tabA = NULL;
152 AjPFeattable tabB = NULL;
153 AjPFeattable outtab = NULL;
154
155 AjPStr seqname = NULL;
156
157 ajint begin;
158 ajint end;
159 float tf;
160
161
162 embInit("twofeat", argc, argv);
163
164 /* input */
165 seqall = ajAcdGetSeqall("sequence");
166
167 /* feature a */
168 asource = ajAcdGetString("asource");
169 atype = ajAcdGetString("atype");
170 asense = ajAcdGetListSingle("asense");
171
172 tf = ajAcdGetFloat("aminscore");
173 aminscore = (ajint) tf;
174 tf = ajAcdGetFloat("amaxscore");
175 amaxscore = (ajint) tf;
176
177 atag = ajAcdGetString("atag");
178 avalue = ajAcdGetString("avalue");
179
180 /* feature b */
181 bsource = ajAcdGetString("bsource");
182 btype = ajAcdGetString("btype");
183 bsense = ajAcdGetListSingle("bsense");
184
185 tf = ajAcdGetFloat("bminscore");
186 bminscore = (ajint) tf;
187 tf = ajAcdGetFloat("bmaxscore");
188 bmaxscore = (ajint) tf;
189
190 btag = ajAcdGetString("btag");
191 bvalue = ajAcdGetString("bvalue");
192
193 /* relation */
194 overlap = ajAcdGetListSingle("overlap");
195 minrange = ajAcdGetInt("minrange");
196 maxrange = ajAcdGetInt("maxrange");
197 rangetype = ajAcdGetListSingle("rangetype");
198 sense = ajAcdGetListSingle("sense");
199 order = ajAcdGetListSingle("order");
200
201 /* output */
202 twoout = ajAcdGetToggle("twoout");
203 typeout = ajAcdGetString("typeout");
204 report = ajAcdGetReport("outfile");
205
206 /* convert feature sense to integer */
207 if(ajStrMatchC(asense, "+"))
208 asensei = +1;
209 else if(ajStrMatchC(asense, "+"))
210 asensei = -1;
211 else
212 asensei = 0;
213
214
215 if(ajStrMatchC(bsense, "+"))
216 bsensei = +1;
217 else if(ajStrMatchC(bsense, "+"))
218 bsensei = -1;
219 else
220 bsensei = 0;
221
222
223 /* convert relation codes to integer values */
224 overlapi = twofeat_get_overlap_type(overlap);
225 rangetypei = twofeat_get_range_type(rangetype);
226 sensei = twofeat_get_sense_type(sense);
227 /* orderi = twofeat_get_order_type(order); Unused */
228
229 seqname = ajStrNew();
230
231 while(ajSeqallNext(seqall, &seq))
232 {
233
234 ajStrAssignC(&seqname, ajSeqGetNameC(seq));
235 begin = ajSeqallGetseqBegin(seqall);
236 end = ajSeqallGetseqEnd(seqall);
237
238 /* make new feature table for A */
239 if(!tabA)
240 tabA = ajFeattableNewSeq(seq);
241
242 /* make new feature table for B */
243 if(!tabB)
244 tabB = ajFeattableNewSeq(seq);
245
246
247 /* go through seq's features adding those that match A to table A */
248 twofeat_find_features(seq, tabA, begin, end, asource, atype, asensei,
249 (float)aminscore, (float)amaxscore, atag,
250 avalue);
251
252 /* go through seq's features adding those that match B to table B */
253 twofeat_find_features(seq, tabB, begin, end, bsource, btype, bsensei,
254 (float)bminscore, (float)bmaxscore, btag,
255 bvalue);
256
257
258 ajDebug("No of hits in tabA: %d\n", ajFeattableGetSize(tabA));
259 ajDebug("No of hits in tabB: %d\n", ajFeattableGetSize(tabB));
260
261 if(ajFeattableGetSize(tabA) && ajFeattableGetSize(tabB))
262 {
263 /* initialise the output feature table */
264 outtab = ajFeattableNewSeq(seq);
265
266
267 /*
268 ** find pairs of hits within the required distance and output
269 ** the results
270 */
271 twofeat_rippledown(tabA, tabB, overlapi, minrange, maxrange,
272 rangetypei, sensei, outtab, twoout,
273 typeout);
274
275 /* write features and tidy up */
276 ajReportWrite(report, outtab, seq);
277 ajFeattableDel(&outtab);
278 }
279
280
281 ajDebug("ajFeattableDel(&tabA)\n");
282 ajFeattableDel(&tabA);
283 ajDebug("ajFeattableDel(&tabB)\n");
284 ajFeattableDel(&tabB);
285 }
286 ajReportSetSeqstats(report, seqall);
287
288 ajStrDel(&seqname);
289 ajSeqDel(&seq);
290 ajSeqallDel(&seqall);
291
292 ajStrDel(&asource);
293 ajStrDel(&atype);
294 ajStrDel(&asense);
295 ajStrDel(&atag);
296 ajStrDel(&avalue);
297
298 ajStrDel(&bsource);
299 ajStrDel(&btype);
300 ajStrDel(&bsense);
301 ajStrDel(&btag);
302 ajStrDel(&bvalue);
303
304 ajStrDel(&overlap);
305 ajStrDel(&rangetype);
306 ajStrDel(&sense);
307 ajStrDel(&order);
308 ajStrDel(&typeout);
309
310 ajReportClose(report);
311 ajReportDel(&report);
312
313 embExit();
314
315 return 0;
316 }
317
318
319
320
321 /* @funcstatic twofeat_rippledown ********************************************
322 **
323 ** Go down the two lists of matches looking for hits within the required
324 ** maximum distance.
325 **
326 ** Foreach feature in tabA
327 ** look at each feature in tabB
328 ** if all the match criteria are valid
329 ** add the matching pair to the list of hits
330 **
331 ** Add any required hits from the hit-list to the output feature table
332 **
333 **
334 ** @param [r] tabA [const AjPFeattable] table A of features to compare to tabB
335 ** @param [r] tabB [const AjPFeattable] table B of features to compare to tabA
336 ** @param [r] overlapi [ajint] types of overlap allowed
337 ** @param [r] minrange [ajint] min distance allowed
338 ** @param [r] maxrange [ajint] max distance allowed
339 ** @param [r] rangetypei [ajint] where to measure the distance from
340 ** @param [r] sensei [ajint] sense relationships allowed
341 ** @param [u] outtab [AjPFeattable] output feature table
342 ** @param [r] twoout [AjBool] True=write both features, else make a single one
343 ** @param [r] typeout [const AjPStr] if a single feature, this is its type name
344 ** @return [void]
345 ** @@
346 ******************************************************************************/
347
twofeat_rippledown(const AjPFeattable tabA,const AjPFeattable tabB,ajint overlapi,ajint minrange,ajint maxrange,ajint rangetypei,ajint sensei,AjPFeattable outtab,AjBool twoout,const AjPStr typeout)348 static void twofeat_rippledown(const AjPFeattable tabA,
349 const AjPFeattable tabB,
350 ajint overlapi, ajint minrange,
351 ajint maxrange, ajint rangetypei,
352 ajint sensei,
353 AjPFeattable outtab, AjBool twoout,
354 const AjPStr typeout)
355 {
356
357
358 AjIList iterA = NULL;
359 AjPFeature gfA = NULL;
360
361 AjIList iterB = NULL;
362 AjPFeature gfB = NULL;
363
364 AjPList hitlist;
365 PHit detail = NULL;
366
367
368 hitlist = ajListNew();
369
370 if(ajFeattableGetSize(tabA))
371 {
372 /* For all features in tabA ... */
373 iterA = ajListIterNewread(tabA->Features);
374 while(!ajListIterDone(iterA))
375 {
376 gfA = ajListIterGet(iterA);
377 ajDebug("In rippledown gfA=%S %d..%d\n",
378 ajFeatGetType(gfA), ajFeatGetStart(gfA),
379 ajFeatGetEnd(gfA));
380
381 /* For all features in tabB ... */
382 if(ajFeattableGetSize(tabB))
383 {
384 iterB = ajListIterNewread(tabB->Features);
385 while(!ajListIterDone(iterB))
386 {
387 gfB = ajListIterGet(iterB);
388 ajDebug("In rippledown gfB=%S %d..%d\n",
389 ajFeatGetType(gfB), ajFeatGetStart(gfB),
390 ajFeatGetEnd(gfB));
391
392 /*
393 ** check for overlap, minrange, maxrange, rangetype,
394 ** sense, order
395 */
396 if(twofeat_check_match(gfA, gfB, &detail, overlapi,
397 minrange, maxrange, rangetypei,
398 sensei))
399 /* push details on hitlist */
400 ajListPush(hitlist, detail);
401 }
402 ajListIterDel(&iterB);
403 }
404 }
405 ajListIterDel(&iterA);
406 }
407
408 /* Put hits in outtab */
409 twofeat_report_hits(hitlist, twoout, typeout, outtab);
410
411
412 ajListFree(&hitlist);
413
414 return;
415 }
416
417
418
419
420 /* @funcstatic twofeat_report_hits ***********************************
421 **
422 ** Outputs the pairs of hits to the output feature table
423 **
424 ** @param [r] hitlist [const AjPList] list of matches (PHit)
425 ** @param [r] twoout [AjBool] True if want pairs of features output
426 ** @param [r] typeout [const AjPStr] name of type if want single features
427 ** @param [u] outtab [AjPFeattable] output feature table
428 ** @return [void]
429 ** @@
430 ******************************************************************************/
431
twofeat_report_hits(const AjPList hitlist,AjBool twoout,const AjPStr typeout,AjPFeattable outtab)432 static void twofeat_report_hits(const AjPList hitlist, AjBool twoout,
433 const AjPStr typeout, AjPFeattable outtab)
434 {
435 char strand;
436 ajint frame = 0;
437 float score = 0.0;
438 AjPStr source = NULL;
439 AjPStr type = NULL;
440 AjPFeature feature;
441 static AjPStr tmp = NULL;
442 ajint start;
443 ajint end;
444
445 AjIList iter = NULL;
446 PHit detail = NULL;
447
448 source = ajStrNew();
449 type = ajStrNew();
450
451 ajStrAssignC(&source,"twofeat");
452 ajStrAssignS(&type, typeout);
453
454
455 iter = ajListIterNewread(hitlist);
456 while(!ajListIterDone(iter))
457 {
458 detail = ajListIterGet(iter);
459
460 if(twoout)
461 {
462 ajFeattableAdd(outtab, ajFeatNewFeat(detail->gfA));
463 ajFeattableAdd(outtab, ajFeatNewFeat(detail->gfB));
464
465 }
466 else
467 {
468 start = detail->Start;
469 end = detail->End;
470
471 /* if both features are -ve then output this, else +ve */
472 if((detail->gfA)->Strand == '-' &&
473 (detail->gfB)->Strand == '-')
474 strand = '-';
475 else
476 strand = '+';
477
478 feature = ajFeatNew(outtab, source, type, start, end,
479 score, strand, frame);
480
481 ajFmtPrintS(&tmp, "*startA %d", ajFeatGetStart(detail->gfA));
482 ajFeatTagAddSS(feature, NULL, tmp);
483
484 ajFmtPrintS(&tmp, "*endA %d", ajFeatGetEnd(detail->gfA));
485 ajFeatTagAddSS(feature, NULL, tmp);
486
487 ajFmtPrintS(&tmp, "*startB %d", ajFeatGetStart(detail->gfB));
488 ajFeatTagAddSS(feature, NULL, tmp);
489
490 ajFmtPrintS(&tmp, "*endB %d", ajFeatGetEnd(detail->gfB));
491 ajFeatTagAddSS(feature, NULL, tmp);
492
493 ajStrDel(&tmp);
494 }
495
496 /* delete hit */
497 twofeat_HitsDel(&detail);
498 }
499 ajListIterDel(&iter);
500
501 ajStrDel(&source);
502 ajStrDel(&type);
503
504 return;
505 }
506
507
508
509
510 /* @funcstatic twofeat_find_features ***********************************
511 **
512 ** Finds seq features matching the required criteria and puts them in tab
513 **
514 ** @param [r] seq [const AjPSeq] the sequence
515 ** @param [u] tab [AjPFeattable] Feature table to populate
516 ** @param [r] begin [ajint] start position in sequence
517 ** @param [r] end [ajint] end position in sequence
518 ** @param [r] source [const AjPStr] source criterion
519 ** @param [r] type [const AjPStr] type criterion
520 ** @param [r] sense [ajint] sense criterion
521 ** @param [r] minscore [float] min score criterion
522 ** @param [r] maxscore [float] max score criterion
523 ** @param [r] tag [const AjPStr] tag criterion
524 ** @param [r] value [const AjPStr] tag value criterion
525 ** @return [void]
526 ** @@
527 ******************************************************************************/
528
twofeat_find_features(const AjPSeq seq,AjPFeattable tab,ajint begin,ajint end,const AjPStr source,const AjPStr type,ajint sense,float minscore,float maxscore,const AjPStr tag,const AjPStr value)529 static void twofeat_find_features(const AjPSeq seq, AjPFeattable tab,
530 ajint begin, ajint end, const AjPStr source,
531 const AjPStr type, ajint sense,
532 float minscore, float maxscore,
533 const AjPStr tag, const AjPStr value)
534 {
535
536 /* get feature table of sequence */
537 const AjPFeattable seqtab = ajSeqGetFeat(seq);
538
539 AjIList iter = NULL;
540 AjPFeature gf = NULL;
541 AjPFeature gfcopy = NULL;
542 AjBool tagsmatch;
543
544
545 tagsmatch = ajFalse;
546
547 /* For all features... */
548 if(ajFeattableGetSize(seqtab))
549 {
550 iter = ajListIterNewread(seqtab->Features);
551 while(!ajListIterDone(iter))
552 {
553 gf = ajListIterGet(iter);
554
555 /* is this feature local and in the sequence range? */
556 if(! ajFeatIsLocalRange(gf, begin, end))
557 continue;
558
559 /* do criterion match */
560 if(twofeat_MatchFeature(gf, source, type, sense, minscore,
561 maxscore, tag, value, &tagsmatch))
562 {
563 ajDebug("Found feature source=%S type=%S %d..%d\n",
564 ajFeatGetSource(gf), ajFeatGetType(gf),
565 ajFeatGetStart(gf), ajFeatGetEnd(gf));
566 /*
567 ** could explicitly make a new feature like this, but it
568 ** is probably safer to let ajFeatNewFeat do it automatically
569 ** gfcopy = ajFeatNew(tab, gf->Source, gf->Type, gf->Start,
570 ** gf->End, gf->Score, gf->Strand, gf->Frame);
571 */
572 gfcopy = ajFeatNewFeat(gf);
573 ajFeattableAdd(tab, gfcopy);
574
575 /* ajFeatTrace(gfcopy); */
576 }
577 }
578 ajListIterDel(&iter);
579 }
580
581 ajDebug("(In twofeat_find_features) No of hits in tab: %d\n",
582 ajFeattableGetSize(tab));
583
584 return;
585 }
586
587
588
589
590 /* @funcstatic twofeat_MatchFeature *****************************************
591 **
592 ** Test if a feature matches a set of criteria
593 **
594 ** @param [r] gf [const AjPFeature] Feature to test
595 ** @param [r] source [const AjPStr] Required Source pattern
596 ** @param [r] type [const AjPStr] Required Type pattern
597 ** @param [r] sense [ajint] Required Sense pattern +1,0,-1 (or other value$
598 ** @param [r] minscore [float] Min required Score pattern
599 ** @param [r] maxscore [float] Max required Score pattern
600 ** @param [r] tag [const AjPStr] Required Tag pattern
601 ** @param [r] value [const AjPStr] Required Value pattern
602 ** @param [w] tagsmatch [AjBool*] true if a join has matching tag/values
603 ** @return [AjBool] True if feature matches criteria
604 ** @@
605 ******************************************************************************/
606
twofeat_MatchFeature(const AjPFeature gf,const AjPStr source,const AjPStr type,ajint sense,float minscore,float maxscore,const AjPStr tag,const AjPStr value,AjBool * tagsmatch)607 static AjBool twofeat_MatchFeature(const AjPFeature gf,
608 const AjPStr source, const AjPStr type,
609 ajint sense, float minscore,
610 float maxscore,
611 const AjPStr tag, const AjPStr value,
612 AjBool *tagsmatch)
613 {
614 AjPStrTok tokens = NULL;
615 AjPStr key = NULL;
616 AjBool val = ajFalse;
617 AjBool scoreok;
618
619 scoreok = (minscore < maxscore);
620
621 /*
622 ** is this a child of a join() ?
623 ** if it is a child, then we use the previous result of MatchPatternTags
624 */
625 if(!ajFeatIsMultiple(gf))
626 *tagsmatch = twofeat_MatchPatternTags(gf, tag, value);
627
628 /* ignore remote IDs */
629 if(!ajFeatIsLocal(gf))
630 return ajFalse;
631
632 /* check source, type, sense, score, tags, values
633 ** Match anything:
634 ** for strings, '*'
635 ** for sense, 0
636 ** for score, maxscore <= minscore
637 */
638
639 if(!embMiscMatchPatternDelimC(ajFeatGetSource(gf), source,",;|") ||
640 (ajFeatGetStrand(gf) == '+' && sense == -1) ||
641 (ajFeatGetStrand(gf) == '-' && sense == +1) ||
642 (scoreok && ajFeatGetScore(gf) < minscore) ||
643 (scoreok && ajFeatGetScore(gf) > maxscore) ||
644 !*tagsmatch)
645 return ajFalse;
646
647 if(ajStrGetLen(type))
648 {
649 val = ajFalse;
650 tokens = ajStrTokenNewC(type, " \t\n\r,;|");
651
652 while (ajStrTokenNextParse(tokens, &key))
653 {
654 if (ajFeatTypeMatchWildS(gf, key))
655 {
656 val = ajTrue;
657 break;
658 }
659 }
660
661 ajStrTokenDel( &tokens);
662 ajStrDel(&key);
663 if(!val)
664 return ajFalse;
665 }
666
667 return ajTrue;
668 }
669
670
671
672
673 /* @funcstatic twofeat_MatchPatternTags **************************************
674 **
675 ** Checks for a match of the tagpattern and valuepattern to at least one
676 ** tag=value pair
677 **
678 ** @param [r] feat [const AjPFeature] Feature to process
679 ** @param [r] tpattern [const AjPStr] tags pattern to match with
680 ** @param [r] vpattern [const AjPStr] values pattern to match with
681 **
682 ** @return [AjBool] ajTrue = found a match
683 ** @@
684 ******************************************************************************/
685
twofeat_MatchPatternTags(const AjPFeature feat,const AjPStr tpattern,const AjPStr vpattern)686 static AjBool twofeat_MatchPatternTags(const AjPFeature feat,
687 const AjPStr tpattern,
688 const AjPStr vpattern)
689 {
690 AjIList titer; /* iterator for feat */
691 static AjPStr tagnam = NULL; /* tag structure */
692 static AjPStr tagval = NULL; /* tag structure */
693 AjBool val = ajFalse; /* returned value */
694 AjBool tval; /* tags result */
695 AjBool vval; /* value result */
696
697
698 /*
699 ** if there are no tags to match, but the patterns are
700 ** both '*', then allow this as a match
701 */
702 if(!ajStrCmpC(tpattern, "*") &&
703 !ajStrCmpC(vpattern, "*"))
704 return ajTrue;
705
706 /* iterate through the tags and test for match to patterns */
707 titer = ajFeatTagIter(feat);
708 while(ajFeatTagval(titer, &tagnam, &tagval))
709 {
710 tval = embMiscMatchPatternDelimC(tagnam, tpattern,",;|");
711 /*
712 ** If tag has no value then
713 ** If vpattern is '*' the value pattern is a match
714 ** Else check vpattern
715 */
716
717 if(!ajStrGetLen(tagval))
718 {
719 if(!ajStrCmpC(vpattern, "*"))
720 vval = ajTrue;
721 else
722 vval = ajFalse;
723 }
724 else
725 vval = embMiscMatchPatternDelimC(tagval, vpattern,",;|");
726
727 if(tval && vval)
728 {
729 val = ajTrue;
730 break;
731 }
732 }
733 ajListIterDel(&titer);
734
735 return val;
736 }
737
738
739
740
741 /* @funcstatic twofeat_HitsNew *****************************************
742 **
743 ** Creates a new PHit object
744 **
745 ** @return [PHit] Hit object
746 ** @@
747 ******************************************************************************/
748
twofeat_HitsNew(void)749 static PHit twofeat_HitsNew(void)
750 {
751 PHit pthis;
752 AJNEW0(pthis);
753
754 return pthis;
755 }
756
757
758
759
760 /* @funcstatic twofeat_HitsDel *****************************************
761 **
762 ** Deletes a PHit object
763 **
764 ** @param [d] pthis [PHit *] Pointer to object to be deleted
765 ** @return [void]
766 ** @@
767 ******************************************************************************/
768
twofeat_HitsDel(PHit * pthis)769 static void twofeat_HitsDel(PHit *pthis)
770 {
771 if(! pthis)
772 return;
773
774 if(! *pthis)
775 return;
776
777 AJFREE(*pthis);
778
779 *pthis = NULL;
780
781 return;
782 }
783
784
785
786
787 /* @funcstatic twofeat_check_match *********************************
788 **
789 ** check for overlap, minrange, maxrange, rangetype, sense, order
790 ** put the results in a PHit object
791 **
792 ** @param [u] gfA [AjPFeature] Feature A - stored for adding new tags
793 ** @param [u] gfB [AjPFeature] Feature B - stored for adding new tags
794 ** @param [w] detail [PHit *] Returned details of match
795 ** @param [r] overlapi [ajint] types of overlap allowed
796 ** @param [r] minrange [ajint] min distance allowed
797 ** @param [r] maxrange [ajint] max distance allowed
798 ** @param [r] rangetypei [ajint] where to measure the distance from
799 ** @param [r] sensei [ajint] sense relationships allowed
800 ** @return [AjBool] True if got a match
801 ** @@
802 ******************************************************************************/
803
twofeat_check_match(AjPFeature gfA,AjPFeature gfB,PHit * detail,ajint overlapi,ajint minrange,ajint maxrange,ajint rangetypei,ajint sensei)804 static AjBool twofeat_check_match(AjPFeature gfA, AjPFeature gfB,
805 PHit *detail, ajint overlapi,
806 ajint minrange, ajint maxrange, ajint
807 rangetypei, ajint sensei)
808 {
809 ajint distance = 0;
810 ajint sA;
811 ajint eA;
812 ajint sB;
813 ajint eB; /* start and end positions */
814 ajint ss;
815 ajint ee;
816 ajint se;
817 ajint es; /* distances from start and end positions */
818 ajint tmp1;
819 ajint tmp2;
820
821 sA = ajFeatGetStart(gfA);
822 eA = ajFeatGetEnd(gfA);
823 sB = ajFeatGetStart(gfB);
824 eB = ajFeatGetEnd(gfB);
825
826 *detail = NULL;
827
828 ajDebug("Next gfA=%S %d..%d\n",
829 ajFeatGetType(gfA), ajFeatGetStart(gfA), ajFeatGetEnd(gfA));
830 ajDebug("Next gfB=%S %d..%d\n",
831 ajFeatGetType(gfB), ajFeatGetStart(gfB), ajFeatGetEnd(gfB));
832
833 /* get distances and absolute distances from each end */
834 ss = abs(sA - sB);
835 ee = abs(eA - eB);
836 se = abs(sA - eB);
837 es = abs(eA - sB);
838
839 switch(rangetypei)
840 {
841 case RA_NEAREST:
842 tmp1 = (ss < ee) ? ss : ee;
843 tmp2 = (se < es) ? se : es;
844 distance = (tmp1 < tmp2) ? tmp1 : tmp2;
845 break;
846
847 case RA_LEFT:
848 distance = ss;
849 break;
850
851 case RA_RIGHT:
852 distance = ee;
853 break;
854
855 case RA_FURTHEST:
856 tmp1 = (ss > ee) ? ss : ee;
857 tmp2 = (se > es) ? se : es;
858 distance = (tmp1 > tmp2) ? tmp1 : tmp2;
859 break;
860
861 default:
862 ajFatal("Unknown range type: %d", rangetypei);
863 }
864
865 ajDebug("Distance: %d\n", distance);
866
867 /* check distance criteria */
868 if(minrange < maxrange)
869 {
870 /* ignore distance if these are the same */
871 ajDebug("Distance < %d: %B\n", minrange, distance < minrange);
872 ajDebug("Distance > %d: %B\n", maxrange, distance > maxrange);
873 if(distance < minrange)
874 return ajFalse;
875
876 if(distance > maxrange)
877 return ajFalse;
878 }
879 ajDebug("Distance OK\n");
880
881 /* check overlap criteria */
882 switch(overlapi)
883 {
884 case OV_ANY:
885 break;
886
887 case OV_OVERLAP: /* want sA in B or sB in A */
888 if((sA<sB || sA>eB) && (sB<sA || sB>eA))
889 return ajFalse;
890 break;
891
892 case OV_NOTOVER: /* want eA before B or sA after B */
893 if((sA>=sB && sA<=eB) || (sB>=sA && sB<=eA))
894 return ajFalse;
895 break;
896
897 case OV_NOTIN: /* want sA or eA out of B and sB or eB out of A */
898 if((sA<=sB && eA>=eB) || (sB<=sA && eB>=eA))
899 return ajFalse;
900 break;
901
902 case OV_AIN: /* want A in B */
903 if(sA<sB || eA>eB)
904 return ajFalse;
905 break;
906
907 case OV_BIN: /* want B in A */
908 if(sB<sA || eB>eA)
909 return ajFalse;
910 break;
911
912 default:
913 ajFatal("Unknown overlap type: %d", overlapi);
914 }
915 ajDebug("Overlap OK\n");
916
917 /* check sense criteria */
918 switch(sensei)
919 {
920 case SN_ANY:
921 break;
922
923 case SN_SAME:
924 if(ajFeatGetStrand(gfA) != ajFeatGetStrand(gfB))
925 return ajFalse;
926 break;
927
928 case SN_OPPOSITE:
929 if(ajFeatGetStrand(gfA) == ajFeatGetStrand(gfB)) return ajFalse;
930
931 default:
932 ajFatal("Unknown sense type: %d", sensei);
933 }
934 ajDebug("Sense OK\n");
935
936 /* check order criteria */
937 /* measure the order from the mid-point of each */
938 switch(sensei)
939 {
940 case SN_ANY:
941 break;
942
943 case OR_AB:
944 tmp1 = (eA-sA)/2;
945 tmp2 = (eB-sB)/2;
946 if(tmp1>tmp2)
947 return ajFalse;
948 break;
949
950 case OR_BA:
951 tmp1 = (eA-sA)/2;
952 tmp2 = (eB-sB)/2;
953 if(tmp1<tmp2)
954 return ajFalse;
955 break;
956
957 default:
958 ajFatal("Unknown sense type: %d", sensei);
959 }
960 ajDebug("Order OK\n");
961
962 /*
963 ** reject any match found to the same feature (same position,
964 ** sense and type)
965 */
966 if(ss == 0 &&
967 ee == 0 &&
968 ajFeatGetStrand(gfA) == ajFeatGetStrand(gfB) &&
969 ajStrMatchS(ajFeatGetType(gfA), ajFeatGetType(gfB)))
970 {
971 ajDebug("Found match of feature to itself\n");
972 return ajFalse;
973 }
974
975 /* if we have a hit, make a PHit object */
976 *detail = twofeat_HitsNew();
977 (*detail)->gfA = gfA;
978 (*detail)->gfB = gfB;
979 (*detail)->distance = distance;
980 (*detail)->Start = (sA < sB) ? sA : sB;
981 (*detail)->End = (eA > eB) ? eA : eB;
982
983 ajDebug("Hit found\n");
984
985 return ajTrue;
986
987 }
988
989
990
991
992 /* @funcstatic twofeat_get_overlap_type *********************************
993 **
994 ** converts the overlap code to an integer
995 **
996 **
997 ** @param [r] overlap [const AjPStr] Overlap code
998 ** @return [ajint] integer value
999 ** @@
1000 ******************************************************************************/
1001
twofeat_get_overlap_type(const AjPStr overlap)1002 static ajint twofeat_get_overlap_type(const AjPStr overlap)
1003 {
1004
1005 if(ajStrMatchC(overlap, "A"))
1006 return OV_ANY;
1007
1008 if(ajStrMatchC(overlap, "O"))
1009 return OV_OVERLAP;
1010
1011 if(ajStrMatchC(overlap, "NO"))
1012 return OV_NOTOVER;
1013
1014 if(ajStrMatchC(overlap, "NW"))
1015 return OV_NOTIN;
1016
1017 if(ajStrMatchC(overlap, "AW"))
1018 return OV_AIN;
1019
1020 if(ajStrMatchC(overlap, "BW")) return OV_BIN;
1021
1022 ajFatal("Unknown Overlap code: %S", overlap);
1023 return -1;
1024 }
1025
1026
1027
1028
1029 /* @funcstatic twofeat_get_range_type *********************************
1030 **
1031 ** converts the range code to an integer
1032 **
1033 **
1034 ** @param [r] rangetype [const AjPStr] Range code
1035 ** @return [ajint] integer value
1036 ** @@
1037 ******************************************************************************/
1038
twofeat_get_range_type(const AjPStr rangetype)1039 static ajint twofeat_get_range_type(const AjPStr rangetype)
1040 {
1041 if(ajStrMatchC(rangetype, "N"))
1042 return RA_NEAREST;
1043
1044 if(ajStrMatchC(rangetype, "L"))
1045 return RA_LEFT;
1046
1047 if(ajStrMatchC(rangetype, "R"))
1048 return RA_RIGHT;
1049
1050 if(ajStrMatchC(rangetype, "F"))
1051 return RA_FURTHEST;
1052
1053 ajFatal("Unknown rangetype code: %S", rangetype);
1054
1055 return -1;
1056 }
1057
1058
1059
1060
1061 /* @funcstatic twofeat_get_sense_type *********************************
1062 **
1063 ** converts the sense code to an integer
1064 **
1065 **
1066 ** @param [r] sense [const AjPStr] sense code
1067 ** @return [ajint] integer value
1068 ** @@
1069 ******************************************************************************/
1070
twofeat_get_sense_type(const AjPStr sense)1071 static ajint twofeat_get_sense_type(const AjPStr sense)
1072 {
1073 if(ajStrMatchC(sense, "A"))
1074 return SN_ANY;
1075
1076 if(ajStrMatchC(sense, "S"))
1077 return SN_SAME;
1078
1079 if(ajStrMatchC(sense, "O"))
1080 return SN_OPPOSITE;
1081
1082 ajFatal("Unknown sense code: %S", sense);
1083
1084 return -1;
1085 }
1086
1087
1088
1089
1090 #if 0
1091 /* @funcstatic twofeat_get_order_type *********************************
1092 **
1093 ** converts the order code to an integer
1094 **
1095 **
1096 ** @param [r] order [const AjPStr] order code
1097 ** @return [ajint] integer value
1098 ** @@
1099 ******************************************************************************/
1100
1101 static ajint twofeat_get_order_type(const AjPStr order)
1102 {
1103 if(ajStrMatchC(order, "A"))
1104 return OR_ANY;
1105
1106 if(ajStrMatchC(order, "AB"))
1107 return OR_AB;
1108
1109 if(ajStrMatchC(order, "BA"))
1110 return OR_BA;
1111
1112 ajFatal("Unknown order code: %S", order);
1113
1114 return -1;
1115 }
1116 #endif
1117