1
2%{
3#include "dnaalign.h"
4#include "hitlist.h"
5
6%}
7
8
9struct DnaMatchPara
10DPRunImpl * dpri
11DnaMatrix * mat
12DnaStartEnd * dse
13int gap
14int ext
15
16%{
17#include "dnamatcher.h"
18
19
20void show_help_DnaMatchPara(FILE * ofp)
21{
22  fprintf(ofp,"Dna Matching Parameters\n");
23  fprintf(ofp,"   -dm_match    [5]   match score\n");
24  fprintf(ofp,"   -dm_mismatch [-4]  mismatch score\n");
25  fprintf(ofp,"   -dm_gappen   [5]   gap open penalty\n");
26  fprintf(ofp,"   -dm_extpen   [1]   gap extension penalty\n");
27
28  show_help_DPRunImpl(ofp);
29}
30
31DnaMatchPara * new_DnaMatchPara_from_argv(int * argc,char ** argv)
32{
33  DnaMatchPara * out;
34  int match    = 5;
35  int mismatch = -10;
36
37  out = DnaMatchPara_alloc();
38
39  strip_out_integer_argument(argc,argv,"dm_match",&match);
40  strip_out_integer_argument(argc,argv,"dm_mismatch",&mismatch);
41
42  assert(mismatch < 0 );
43
44  out->dpri = new_DPRunImpl_from_argv(argc,argv);
45  out->mat = identity_DnaMatrix(match,mismatch);
46  out->dse = DnaStartEnd_from_policy("local");
47
48  out->gap = 30;
49  out->ext = 20;
50
51  strip_out_integer_argument(argc,argv,"dm_gappen",&out->gap);
52  strip_out_integer_argument(argc,argv,"dm_extpen",&out->ext);
53
54
55  return out;
56}
57
58
59HitList * HitList_from_Sequence_SequenceSet_DNA(Sequence * query,SequenceSet * set,DnaMatchPara * p)
60{
61  int i;
62  HitList * out;
63  HitPair * pair;
64  HitAln * aln;
65
66  AlnBlock * forward;
67  AlnBlock * reverse;
68
69  Sequence * rev;
70
71  char buffer[512];
72
73  out = HitList_alloc_std();
74
75  for(i=0;i<set->len;i++) {
76
77    rev = reverse_complement_Sequence(set->set[i]);
78
79    ckfree(rev->name);
80    sprintf(buffer,"%s.reverse",set->set[i]->name);
81    rev->name  = stringalloc(buffer);
82
83    pair = HitPair_alloc_std();
84
85    aln = HitAln_alloc();
86
87    forward = make_align_dnaalign(query,set->set[i],p->mat,p->dse,-p->gap,-p->ext,-p->gap,-p->ext,p->dpri);
88    reverse = make_align_dnaalign(query,rev,p->mat,p->dse,-p->gap,-p->ext,-p->gap,-p->ext,p->dpri);
89
90    if( forward->score > reverse->score ) {
91      pair->query  = hard_link_Sequence(query);
92      pair->target = hard_link_Sequence(set->set[i]);
93      aln->alb = hard_link_AlnBlock(forward);
94    } else {
95      pair->query  = hard_link_Sequence(query);
96      pair->target = hard_link_Sequence(rev);
97      aln->alb = hard_link_AlnBlock(reverse);
98    }
99
100
101
102    add_HitPair(pair,aln);
103
104    aln->raw_score = pair->raw_score = aln->alb->score;
105    add_HitList(out,pair);
106
107
108    free_AlnBlock(forward);
109    free_AlnBlock(reverse);
110
111    free_Sequence(rev);
112
113  }
114
115  return out;
116}
117