ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/seqalign/seqalign.cpp
(Generate patch)
# Line 52 | Line 52
52   return false;
53   }
54  
55
56 GXAlnInfo* match_RightEnd(const char* seqa, int seqa_len, const char* seqb, int seqb_len, SGreedyAlignMem* gxmem) {
57  bool editscript=false;
58  #ifdef GDEBUG
59   editscript=true;
60   GMessage("==========> matching Right (3') end ========\n");
61  #endif
62  GList<GXSeed> rseeds(true,true,false);
63  GXBandSet* alnbands=collectSeeds_R(rseeds, seqa, seqa_len, seqb, seqb_len);
64  GList<GXAlnInfo> galns(true,true,false);
65  int max_top_bands=5;
66  int top_band_count=0;
67  GList<GXSeed> anchor_seeds(cmpSeedDiag, NULL, true); //stores unique seeds per diagonal
68  for (int b=0;b<alnbands->Count();b++) {
69     if (alnbands->Get(b)->score<4) break;
70     top_band_count++;
71     GXBand& band=*(alnbands->Get(b));
72     band.seeds.setSorted(cmpSeedScore);
73     anchor_seeds.Add(band.seeds.First());
74     if (anchor_seeds.Count()>2 || top_band_count>max_top_bands) break;
75     }
76  #ifdef GDEBUG
77  GMessage("::: Collected %d anchor seeds.\n",anchor_seeds.Count());
78  #endif
79
80  for (int i=0;i<anchor_seeds.Count();i++) {
81    GXSeed& aseed=*anchor_seeds[i];
82    int a1=aseed.a_ofs+(aseed.len>>1)+1;
83    int a2=aseed.b_ofs+(aseed.len>>1)+1;
84    GXAlnInfo* alninfo=GreedyAlignRegion(seqa, a1, seqa_len,
85                            seqb, a2, seqb_len, galn_TrimRight, editscript, gxmem,
86                            match_reward, mismatch_penalty, Xdrop);
87    if (alninfo && alninfo->pid>73 && alninfo->sr>seqb_len-4 && alninfo->ql<4)
88             galns.AddIfNew(alninfo, true);
89        else delete alninfo;
90    }
91
92  #ifdef GDEBUG
93  for (int i=0;i<galns.Count();i++) {
94    GXAlnInfo* alninfo=galns[i];
95    GMessage("a(%d..%d) align to b(%d..%d), score=%d, pid=%4.2f\n", alninfo->ql, alninfo->qr,
96                         alninfo->sl, alninfo->sr, alninfo->score, alninfo->pid);
97    if (alninfo->gapinfo!=NULL) {
98      GMessage("Alignment:\n");
99      alninfo->gapinfo->printAlignment(stderr, seqa, seqa_len, seqb, seqb_len);
100      }
101    }
102  #endif
103  //---- done
104  delete alnbands;
105  if (galns.Count()) {
106    GXAlnInfo* bestaln=galns.Shift();
107    return bestaln;
108    }
109  else return NULL;
110 }
111
112 GXAlnInfo* match_LeftEnd(const char* seqa, int seqa_len, const char* seqb, int seqb_len, SGreedyAlignMem* gxmem) {
113  bool editscript=false;
114  #ifdef GDEBUG
115   editscript=true;
116   GMessage("==========> matching Left (5') end ========\n");
117  #endif
118  GList<GXSeed> rseeds(true,true,false);
119  GXBandSet* alnbands=collectSeeds_L(rseeds, seqa, seqa_len, seqb, seqb_len);
120  GList<GXAlnInfo> galns(true,true,false);
121  int max_top_bands=5;
122  int top_band_count=0;
123  GList<GXSeed> anchor_seeds(cmpSeedDiag, NULL, true); //stores unique seeds per diagonal
124  for (int b=0;b<alnbands->Count();b++) {
125     if (alnbands->Get(b)->score<4) break;
126     top_band_count++;
127     GXBand& band=*(alnbands->Get(b));
128     band.seeds.setSorted(cmpSeedScore);
129     anchor_seeds.Add(band.seeds.First());
130     if (anchor_seeds.Count()>2 || top_band_count>max_top_bands) break;
131     }
132  #ifdef GDEBUG
133  GMessage("::: Collected %d anchor seeds.\n",anchor_seeds.Count());
134  #endif
135
136  for (int i=0;i<anchor_seeds.Count();i++) {
137    GXSeed& aseed=*anchor_seeds[i];
138    int a1=aseed.a_ofs+(aseed.len>>1)+1;
139    int a2=aseed.b_ofs+(aseed.len>>1)+1;
140    GXAlnInfo* alninfo=GreedyAlignRegion(seqa, a1, seqa_len,
141                            seqb, a2, seqb_len, galn_TrimLeft, editscript, gxmem,
142                            match_reward, mismatch_penalty, Xdrop-2);
143    if (alninfo && alninfo->pid>82 && alninfo->sl<3 && alninfo->qr>seqa_len-4)
144         //TODO: at 5' end we should use a higher pid threshold
145         // we could also decrease the X-drop value!
146            galns.AddIfNew(alninfo, true);
147       else delete alninfo;
148    }
149
150  #ifdef GDEBUG
151  for (int i=0;i<galns.Count();i++) {
152    GXAlnInfo* alninfo=galns[i];
153    GMessage("a(%d..%d) align to b(%d..%d), score=%d, pid=%4.2f\n", alninfo->ql, alninfo->qr,
154                         alninfo->sl, alninfo->sr, alninfo->score, alninfo->pid);
155    if (alninfo->gapinfo!=NULL) {
156      GMessage("Alignment:\n");
157      alninfo->gapinfo->printAlignment(stderr, seqa, seqa_len, seqb, seqb_len);
158      }
159    }
160  #endif
161  //---- done
162  delete alnbands;
163  if (galns.Count()) {
164    GXAlnInfo* bestaln=galns.Shift();
165    return bestaln;
166    }
167  else return NULL;
168 }
169
170
171
55   void findTrimming(char* aseq, int alen, char* rseq, int rlen, SGreedyAlignMem* gxmem=NULL) {
56  
57 < GXAlnInfo* r_bestaln=match_RightEnd(aseq, alen, rseq, rlen, gxmem);
57 > GXAlnInfo* r_bestaln=match_RightEnd(aseq, alen, rseq, rlen, gxmem, match_reward, mismatch_penalty, Xdrop);
58  
59   #ifdef GDEBUG
60    if (r_bestaln)
# Line 189 | Line 72
72    reverseChars(aseq);
73    reverseChars(rseq);
74  
75 <  GXAlnInfo* l_bestaln=match_LeftEnd(aseq, alen, rseq, rlen, gxmem);
75 >  GXAlnInfo* l_bestaln=match_LeftEnd(aseq, alen, rseq, rlen, gxmem, match_reward, mismatch_penalty, Xdrop);
76  
77   #ifdef GDEBUG
78    if (l_bestaln)

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines