1 /* @source matcher application
2 **
3 ** MATCHER finds the best local alignments between two sequences
4 ** version 2.0u4 Feb. 1996
5 ** Please cite:
6 ** X. Huang and W. Miller (1991) Adv. Appl. Math. 12:373-381
7 **
8 ** @author Copyright (C) Ian Longden (il@sanger.ac.uk)
9 ** @@
10 **
11 ** This program is free software; you can redistribute it and/or
12 ** modify it under the terms of the GNU General Public License
13 ** as published by the Free Software Foundation; either version 2
14 ** of the License, or (at your option) any later version.
15 **
16 ** This program is distributed in the hope that it will be useful,
17 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
18 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 ** GNU General Public License for more details.
20 **
21 ** You should have received a copy of the GNU General Public License
22 ** along with this program; if not, write to the Free Software
23 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
24 ******************************************************************************/
25
26 #include "emboss.h"
27
28
29
30
31 /* @macro matchergap **********************************************************
32 **
33 ** sets k-symbol indel score
34 **
35 ** @param [r] k [ajint] Symbol
36 ** @return [void]
37 ******************************************************************************/
38
39 #define matchergap(k) ((k) <= 0 ? 0 : q+r*(k)) /* k-symbol indel score */
40
41
42
43
44 static ajint tt;
45
46 static ajint *sapp; /* Current script append ptr */
47 static ajint last; /* Last script op appended */
48
49 static ajint I, J; /* current positions of A ,B */
50 static ajint al_len; /* length of alignment */
51
52
53
54
55 /* @macro MATCHERDEL **********************************************************
56 **
57 ** Macro for a "Delete k" operation
58 **
59 ** @param [r] k [ajint] Undocumented
60 ** @return [void]
61 ******************************************************************************/
62
63 #define MATCHERDEL(k) \
64 { I += k; \
65 al_len += k; \
66 if (last < 0) \
67 last = sapp[-1] -= (k); \
68 else \
69 last = *sapp++ = -(k); \
70 }
71
72
73
74
75 /* @macro MATCHERINS **********************************************************
76 **
77 ** Macro for an "Insert k" operation
78 **
79 ** @param [r] k [ajint] Undocumented
80 ** @return [void]
81 ******************************************************************************/
82
83 #define MATCHERINS(k) \
84 { J += k; \
85 al_len += k; \
86 if (last < 0) \
87 { sapp[-1] = (k); *sapp++ = last; } \
88 else \
89 last = *sapp++ = (k); \
90 }
91
92
93
94
95 /* @macro MATCHERREP **********************************************************
96 **
97 ** Macro for a "Replace" operation
98 **
99 ** @return [void]
100 ******************************************************************************/
101
102 #define MATCHERREP \
103 { last = *sapp++ = 0; \
104 al_len += 1; \
105 }
106
107
108
109
110 /* @macro MATCHERDIAG *********************************************************
111 **
112 ** assigns value to x if (ii,jj) is never used before
113 **
114 ** @param [r] ii [ajint] Undocumented
115 ** @param [r] jj [ajint] Undocumented
116 ** @param [r] x [ajint] Undocumented
117 ** @param [r] value [ajint] Undocumented
118 ** @return [void]
119 ******************************************************************************/
120
121 #define MATCHERDIAG(ii, jj, x, value) \
122 { \
123 for (tt = 1, z = row[(ii)]; z != PAIRNULL; z = z->NEXT) \
124 if (z->COL == (jj)) \
125 { tt = 0; break; } \
126 if (tt) \
127 x = (value); \
128 }
129
130
131
132
133 /* @macro MATCHERORDER ********************************************************
134 **
135 ** replace (ss1, xx1, yy1) by (ss2, xx2, yy2) if the latter is large
136 **
137 ** @param [r] ss1 [ajint] Undocumented
138 ** @param [r] xx1 [ajint] Undocumented
139 ** @param [r] yy1 [ajint] Undocumented
140 ** @param [r] ss2 [ajint] Undocumented
141 ** @param [r] xx2 [ajint] Undocumented
142 ** @param [r] yy2 [ajint] Undocumented
143 ** @return [void]
144 ******************************************************************************/
145
146 #define MATCHERORDER(ss1, xx1, yy1, ss2, xx2, yy2) \
147 { if (ss1 < ss2) \
148 { ss1 = ss2; xx1 = xx2; yy1 = yy2; } \
149 else \
150 if (ss1 == ss2) \
151 { if (xx1 < xx2) \
152 { xx1 = xx2; yy1 = yy2; } \
153 else \
154 if (xx1 == xx2 && yy1 < yy2) \
155 yy1 = yy2; \
156 } \
157 }
158
159
160
161
162 static AjPSeq seq;
163 static AjPSeq seq2;
164 static ajint **sub;
165
166
167
168
169 /* @datastatic vertex *********************************************************
170 **
171 ** Matcher internals
172 **
173 ** @alias NODE
174 ** @alias vertexptr
175 **
176 ** @attr SCORE [ajint] Undocumented
177 ** @attr STARI [ajint] Undocumented
178 ** @attr STARJ [ajint] Undocumented
179 ** @attr ENDI [ajint] Undocumented
180 ** @attr ENDJ [ajint] Undocumented
181 ** @attr TOP [ajint] Undocumented
182 ** @attr BOT [ajint] Undocumented
183 ** @attr LEFT [ajint] Undocumented
184 ** @attr RIGHT [ajint] Undocumented
185 ******************************************************************************/
186
187 typedef struct NODE
188 {
189 ajint SCORE;
190 ajint STARI;
191 ajint STARJ;
192 ajint ENDI;
193 ajint ENDJ;
194 ajint TOP;
195 ajint BOT;
196 ajint LEFT;
197 ajint RIGHT;
198 } vertex;
199 #define vertexptr vertex*
200
201
202
203
204 vertexptr *LIST; /* an array for saving k best scores */
205 vertexptr low = 0; /* lowest score node in LIST */
206 vertexptr most = 0; /* latestly accessed node in LIST */
207 static ajint *CC, *DD; /* saving matrix scores */
208 static ajint *RR, *SS, *EE, *FF; /* saving start-points */
209 static ajint *HH, *WW; /* saving matrix scores */
210 static ajint *II, *JJ, *XX, *YY; /* saving start-points */
211 static ajint m1, mm, n1, nn; /* boundaries of recomputed area */
212 static ajint rl, cl; /* left and top boundaries */
213 static ajint lmin; /* minimum score in LIST */
214 static ajint flag; /* indicate if recomputation needed */
215
216 static ajint q, r; /* gap penalties */
217 static ajint qr; /* qr = q + r */
218
219
220
221
222 /* @datastatic pair ***********************************************************
223 **
224 ** Matcher internals
225 **
226 ** @alias ONE
227 ** @alias pairptr
228 **
229 ** @attr NEXT [struct ONE*] Undocumented
230 ** @attr COL [ajint] Undocumented
231 ** @attr Padding [char[4]] Padding to alignment boundary
232 ******************************************************************************/
233
234 typedef struct ONE
235 {
236 struct ONE *NEXT;
237 ajint COL;
238 char Padding[4];
239 } pair;
240 #define pairptr pair*
241
242 pairptr *row; /* for saving used aligned pairs */
243 pairptr z; /* for saving used aligned pairs */
244
245
246
247 #define PAIRNULL (pairptr)NULL
248
249 AjPMatrix matrix = NULL;
250 AjPSeqCvt cvt = NULL;
251
252
253
254
255 static void matcher_Sim(AjPAlign align,
256 const char A[], const char B[],
257 const AjPSeq seq0, const AjPSeq seq1, ajuint K,
258 ajint Q, ajint R, ajint beg, ajint beg2, ajint nseq);
259 static ajint matcher_BigPass(const char A[], const char B[],
260 ajint M, ajint N, ajint K,
261 ajint nseq, ajint *numnode);
262 static ajint matcher_Locate(const char A[], const char B[], ajint nseq,
263 ajint numnode);
264 static ajint matcher_SmallPass(const char A[], const char B[],
265 ajint count, ajint nseq, ajint *numnode);
266 static ajint matcher_Diff(const char A[], const char B[],
267 ajint M, ajint N, ajint tb,
268 ajint te);
269 static ajint matcher_Calcons(char *seqc0, char *seqc1,
270 const ajint *res,
271 ajint min0, ajint min1,
272 ajint max0, ajint max1,
273 ajint *nc, ajint *nident);
274 static ajint matcher_Addnode(ajint c, ajint ci, ajint cj, ajint i, ajint j,
275 ajint K, ajint cost, ajint *numnode);
276 static ajint matcher_NoCross(ajint numnode);
277 static vertexptr matcher_Findmax(ajint *numnode);
278
279
280
281
282 /* @prog matcher **************************************************************
283 **
284 ** Finds the best local alignments between two sequences
285 **
286 ******************************************************************************/
287
main(int argc,char ** argv)288 int main(int argc, char **argv)
289 {
290 AjPStr aa0str = 0;
291 AjPStr aa1str = 0;
292 const char *s1;
293 const char *s2;
294 ajint gdelval;
295 ajint ggapval;
296 ajuint i;
297 ajint K;
298
299 AjPAlign align = NULL;
300
301 embInit("matcher", argc, argv);
302
303 seq = ajAcdGetSeq("asequence");
304 ajSeqTrim(seq);
305 seq2 = ajAcdGetSeq("bsequence");
306 ajSeqTrim(seq2);
307 matrix = ajAcdGetMatrix("datafile");
308 K = ajAcdGetInt("alternatives");
309 gdelval = ajAcdGetInt("gapopen");
310 ggapval = ajAcdGetInt("gapextend");
311 align = ajAcdGetAlign("outfile");
312
313
314 /*
315 create sequence indices. i.e. A->0, B->1 ... Z->25 etc.
316 This is done so that ajAZToInt has only to be done once for
317 each residue in the sequence
318 */
319
320 ajSeqFmtUpper(seq);
321 ajSeqFmtUpper(seq2);
322
323 s1 = ajStrGetPtr(ajSeqGetSeqS(seq));
324 s2 = ajStrGetPtr(ajSeqGetSeqS(seq2));
325
326 sub = ajMatrixGetMatrix(matrix);
327 cvt = ajMatrixGetCvt(matrix);
328
329
330 aa0str = ajStrNewRes(2+ajSeqGetLen(seq)); /* length + blank + trailing null */
331 aa1str = ajStrNewRes(2+ajSeqGetLen(seq2));
332 ajStrAppendK(&aa0str,' ');
333 ajStrAppendK(&aa1str,' ');
334
335 for(i=0;i<ajSeqGetLen(seq);i++)
336 ajStrAppendK(&aa0str,(char)ajSeqcvtGetCodeK(cvt, *s1++));
337
338 for(i=0;i<ajSeqGetLen(seq2);i++)
339 ajStrAppendK(&aa1str,ajSeqcvtGetCodeK(cvt, *s2++));
340
341 matcher_Sim(align, ajStrGetPtr(aa0str),ajStrGetPtr(aa1str),
342 seq,seq2,
343 K,(gdelval-ggapval),ggapval,
344 ajSeqGetOffset(seq), ajSeqGetOffset(seq2), 2);
345
346 ajStrDel(&aa0str);
347 ajStrDel(&aa1str);
348
349 ajSeqDel(&seq);
350 ajSeqDel(&seq2);
351
352 embExit();
353
354 return 0;
355 }
356
357
358
359
360 /* @funcstatic matcher_Sim ****************************************************
361 **
362 ** Smith Waterman Eggert local alignment
363 **
364 ** @param [u] align [AjPAlign] Alignment object
365 ** @param [r] A [const char*] Sequence A with trailing blank
366 ** @param [r] B [const char*] Sequence B with trailing blank
367 ** @param [r] seq0 [const AjPSeq] Sequence A
368 ** @param [r] seq1 [const AjPSeq] Sequence B
369 ** @param [r] K [ajuint] Number of best alignments to report
370 ** @param [r] Q [ajint] Gap penalty (minus extension penalty)
371 ** @param [r] R [ajint] Gap extension penalty
372 ** @param [r] beg0 [ajint] Offset of sequence A
373 ** @param [r] beg1 [ajint] Offset of sequence B
374 ** @param [r] nseq [ajint] Number of sequences
375 ** @return [void]
376 ******************************************************************************/
377
matcher_Sim(AjPAlign align,const char * A,const char * B,const AjPSeq seq0,const AjPSeq seq1,ajuint K,ajint Q,ajint R,ajint beg0,ajint beg1,ajint nseq)378 static void matcher_Sim(AjPAlign align,
379 const char *A,const char *B,
380 const AjPSeq seq0, const AjPSeq seq1,ajuint K,ajint Q,
381 ajint R, ajint beg0, ajint beg1, ajint nseq)
382 {
383 ajint endi;
384 ajint endj;
385 ajint stari;
386 ajint starj; /* endpoint and startpoint */
387 ajint score; /* the max score in LIST */
388 ajint count; /* maximum size of list */
389 register ajuint i;
390 register ajuint j; /* row and column indices */
391 ajint *S; /* saving operations for diff */
392 ajint nc;
393 ajint nident; /* for display */
394 vertexptr cur; /* temporary pointer */
395 ajint seq0len;
396 ajint seq1len;
397
398 AjPSeq res0 = NULL;
399 AjPSeq res1 = NULL;
400
401 char *seqc0, *seqc1; /* aligned sequences */
402
403 ajint numnode; /* the number of nodes in LIST */
404
405 ajint min0;
406 ajint min1;
407 ajint max0;
408 ajint max1;
409
410 seq0len = ajSeqGetLen(seq0);
411 seq1len = ajSeqGetLen(seq1);
412
413 /* allocate space for consensus */
414 i = (seq0len+seq1len+1);
415 AJCNEW(seqc0,i);
416 AJCNEW(seqc1,i);
417 /* allocate space for all vectors */
418
419 j = (seq1len + 1) /* * sizeof(ajint)*/;
420 AJCNEW(CC, j);
421 AJCNEW(DD, j);
422 AJCNEW(RR, j);
423 AJCNEW(SS, j);
424 AJCNEW(EE, j);
425 AJCNEW(FF, j);
426
427 i = (seq0len + 1) /* * sizeof(ajint)*/;
428 AJCNEW(HH, i);
429 AJCNEW(WW, i);
430 AJCNEW(II, i);
431 AJCNEW(JJ, i);
432 AJCNEW(XX, i);
433 AJCNEW(YY, i);
434 AJCNEW(S, (i+j));
435 AJCNEW0(row,(seq0len + 1));
436
437 /* set up list for each row (already zeroed by AJCNEW0 macro) */
438 /* for(i = 1; i <= M; i++) row[i]= PAIRNULL;*/
439
440 /* vv = *sub[0];*/
441 q = Q;
442 r = R;
443 qr = q + r;
444
445 AJCNEW(LIST, K);
446 for(i = 0; i < K; i++)
447 AJNEW0(LIST[i]);
448
449 numnode = lmin = 0;
450 matcher_BigPass(A,B,seq0len,seq1len,K,nseq, &numnode);
451
452 /* Report the K best alignments one by one. After each alignment is
453 output, recompute part of the matrix. First determine the size
454 of the area to be recomputed, then do the recomputation */
455
456 for(count = K - 1; count >= 0; count--)
457 {
458 if(numnode == 0)
459 ajFatal("The number of alignments computed is too large");
460 cur = matcher_Findmax(&numnode);
461 score = cur->SCORE;
462 stari = ++cur->STARI;
463 starj = ++cur->STARJ;
464 endi = cur->ENDI;
465 endj = cur->ENDJ;
466 m1 = cur->TOP;
467 mm = cur->BOT;
468 n1 = cur->LEFT;
469 nn = cur->RIGHT;
470 rl = endi - stari + 1;
471 cl = endj - starj + 1;
472 I = stari - 1;
473 J = starj - 1;
474 sapp = S;
475 last = 0;
476 al_len = 0;
477 matcher_Diff(&A[stari]-1, &B[starj]-1,rl,cl,q,q);
478
479 min0 = stari-1;
480 min1 = starj-1;
481 max0 = stari+rl-1;
482 max1 = starj+cl-1;
483 matcher_Calcons(seqc0,seqc1,
484 S, min0, min1, max0, max1, &nc, &nident);
485 /* percent = (double)nident*100.0/(double)nc; Unused */
486
487 ajDebug("Matcher min: %d %d max: %d %d beg: %d %d nc: %d\n",
488 min0, min1, max0, max1, beg0, beg1, nc);
489 ajDebug("Matcher offsets: %d %d end: %d %d len: %d %d alen: %d %d\n",
490 ajSeqGetOffset(seq0), ajSeqGetOffset(seq1),
491 ajSeqGetOffend(seq0), ajSeqGetOffend(seq1),
492 ajSeqGetLen(seq0), ajSeqGetLen(seq1),
493 strlen(seqc0), strlen(seqc1));
494 ajDebug("Matcher seqc0: %s\n", seqc0);
495 ajDebug("Matcher seqc1: %s\n", seqc1);
496
497 res0 = ajSeqNewRangeC(seqc0, min0+beg0, ajSeqGetOffend(seq0) + ajSeqGetLen(seq0) - max0,
498 ajSeqIsReversed(seq0));
499 ajSeqAssignUsaS(res0, ajSeqGetUsaS(seq0));
500 ajSeqAssignNameS(res0, ajSeqGetNameS(seq0));
501
502 res1 = ajSeqNewRangeC(seqc1, min1+beg1, ajSeqGetOffend(seq1) + ajSeqGetLen(seq1) - max1,
503 ajSeqIsReversed(seq1));
504 ajSeqAssignNameS(res1, ajSeqGetNameS(seq1));
505 ajSeqAssignUsaS(res1, ajSeqGetUsaS(seq1));
506
507 ajAlignDefineSS(align, res0, res1);
508 ajSeqDel(&res0);
509 ajSeqDel(&res1);
510
511 ajAlignSetGapI(align, Q+R, R);
512 ajAlignSetScoreI(align, score);
513 ajAlignSetMatrixInt(align, matrix);
514 ajAlignSetStats(align, -1, nc, nident, -1, -1, NULL);
515
516 ajAlignTrace(align);
517
518 if(count)
519 {
520 flag = 0;
521 matcher_Locate(A,B,nseq, numnode);
522 if(flag)
523 matcher_SmallPass(A,B,count,nseq, &numnode);
524 }
525 }
526
527 ajAlignWrite(align);
528
529 /* now free all the memory */
530
531 ajAlignClose(align);
532
533 ajAlignDel(&align);
534
535 AJFREE(CC);
536 AJFREE(DD);
537 AJFREE(RR);
538 AJFREE(SS);
539 AJFREE(EE);
540 AJFREE(FF);
541 AJFREE(HH);
542 AJFREE(WW);
543 AJFREE(II);
544 AJFREE(JJ);
545 AJFREE(XX);
546 AJFREE(YY);
547 AJFREE(S);
548
549 for(i=0; i<=ajSeqGetLen(seq);i++)
550 {
551 pairptr this;
552 pairptr next;
553 if(row[i])
554 {
555 this = row[i];
556 next = this->NEXT;
557 while(next)
558 {
559 AJFREE(this);
560 this = next;
561 next= this->NEXT;
562 }
563 AJFREE(this);
564 }
565 }
566
567 AJFREE(row);
568 for(i = 0; i < K; i++)
569 AJFREE(LIST[i]);
570 AJFREE(LIST);
571
572
573 AJFREE(seqc0);
574 AJFREE(seqc1);
575
576 return;
577 }
578
579
580
581
582 /* @funcstatic matcher_NoCross ************************************************
583 **
584 ** return 1 if no node in LIST shares vertices with the area
585 **
586 ** @param [r] numnode [ajint] Undocumented
587 ** @return [ajint] 1 if no node shares vertices
588 ******************************************************************************/
589
matcher_NoCross(ajint numnode)590 static ajint matcher_NoCross(ajint numnode)
591 {
592 vertexptr cur;
593 register ajint i;
594
595 for(i = 0; i < numnode; i++)
596 {
597 cur = LIST[i];
598 if(cur->STARI <= mm && cur->STARJ <= nn && cur->BOT >= m1-1 &&
599 cur->RIGHT >= n1-1 && ( cur->STARI < rl || cur->STARJ < cl))
600 {
601 if(cur->STARI < rl)
602 rl = cur->STARI;
603
604 if(cur->STARJ < cl)
605 cl = cur->STARJ;
606 flag = 1;
607 break;
608 }
609 }
610
611 if(i == numnode)
612 return 1;
613
614 return 0;
615 }
616
617
618
619
620 /* @funcstatic matcher_BigPass ************************************************
621 **
622 ** Undocumented
623 **
624 ** @param [r] A [const char*] Undocumented
625 ** @param [r] B [const char*] Undocumented
626 ** @param [r] M [ajint] Undocumented
627 ** @param [r] N [ajint] Undocumented
628 ** @param [r] K [ajint] Undocumented
629 ** @param [r] nseq [ajint] Number of sequences
630 ** @param [w] numnode [ajint*] Number of alignments in list
631 ** @return [ajint] Undocumented
632 ******************************************************************************/
633
matcher_BigPass(const char * A,const char * B,ajint M,ajint N,ajint K,ajint nseq,ajint * numnode)634 static ajint matcher_BigPass(const char *A,const char *B,
635 ajint M,ajint N,ajint K,
636 ajint nseq, ajint *numnode)
637 {
638 register ajint i;
639 register ajint j; /* row and column indices */
640 register ajint c; /* best score at current point */
641 register ajint f; /* best score ending with insertion */
642 register ajint d; /* best score ending with deletion */
643 register ajint p; /* best score at (i-1, j-1) */
644 register ajint ci;
645 register ajint cj; /* end-point associated with c */
646 register ajint di;
647 register ajint dj; /* end-point associated with d */
648 register ajint fi;
649 register ajint fj; /* end-point associated with f */
650 register ajint pi;
651 register ajint pj; /* end-point associated with p */
652 ajint *va; /* pointer to vv(A[i], B[j]) */
653
654 /*
655 ** Compute the matrix and save the top K best scores in LIST
656 ** CC : the scores of the current row
657 ** RR and EE : the starting point that leads to score CC
658 ** DD : the scores of the current row, ending with deletion
659 ** SS and FF : the starting point that leads to score DD
660 */
661
662 /* Initialize the 0 th row */
663 for(j = 1; j <= N; j++)
664 {
665 CC[j] = 0;
666 RR[j] = 0;
667 EE[j] = j;
668 DD[j] = - (q);
669 SS[j] = 0;
670 FF[j] = j;
671 }
672
673 for(i = 1; i <= M; i++)
674 {
675 c = 0; /* Initialize column 0 */
676 f = - (q);
677 ci = fi = i;
678 va = sub[(ajint)A[i]];
679
680 if(nseq == 2)
681 {
682 p = 0;
683 pi = i - 1;
684 cj = fj = pj = 0;
685 }
686 else
687 {
688 p = CC[i];
689 pi = RR[i];
690 pj = EE[i];
691 cj = fj = i;
692 }
693
694 for(j = (nseq == 2 ? 1 : (i+1)); j <= N; j++)
695 {
696 f = f - r;
697 c = c - qr;
698 MATCHERORDER(f, fi, fj, c, ci, cj)
699 c = CC[j] - qr;
700 ci = RR[j];
701 cj = EE[j];
702 d = DD[j] - r;
703 di = SS[j];
704 dj = FF[j];
705 MATCHERORDER(d, di, dj, c, ci, cj)
706 c = 0;
707 MATCHERDIAG(i, j, c, p+va[(ajint)B[j]]) /* diagonal */
708
709 if(c <= 0)
710 {
711 c = 0;
712 ci = i;
713 cj = j;
714 }
715 else
716 {
717 ci = pi; cj = pj; }
718 MATCHERORDER(c, ci, cj, d, di, dj)
719 MATCHERORDER(c, ci, cj, f, fi, fj)
720 p = CC[j];
721 CC[j] = c;
722 pi = RR[j];
723 pj = EE[j];
724 RR[j] = ci;
725 EE[j] = cj;
726 DD[j] = d;
727 SS[j] = di;
728 FF[j] = dj;
729 if(c > lmin) /* add the score into list */
730 lmin = matcher_Addnode(c, ci, cj, i, j, K, lmin, numnode);
731 }
732 }
733
734 return 0;
735 }
736
737
738
739
740 /* @funcstatic matcher_Locate *************************************************
741 **
742 ** Undocumented
743 **
744 ** @param [r] A [const char*] Undocumented
745 ** @param [r] B [const char*] Undocumented
746 ** @param [r] nseq [ajint] Number of sequences
747 ** @param [r] numnode [ajint] Number of alignments in list
748 ** @return [ajint] Undocumented
749 ******************************************************************************/
750
matcher_Locate(const char * A,const char * B,ajint nseq,ajint numnode)751 static ajint matcher_Locate(const char *A,const char *B,ajint nseq,
752 ajint numnode)
753 {
754 register ajint i;
755 register ajint j; /* row and column indices */
756 register ajint c; /* best score at current point */
757 register ajint f; /* best score ending with insertion */
758 register ajint d; /* best score ending with deletion */
759 register ajint p; /* best score at (i-1, j-1) */
760 register ajint ci;
761 register ajint cj; /* end-point associated with c */
762 register ajint di=0;
763 register ajint dj=0; /* end-point associated with d */
764 register ajint fi;
765 register ajint fj; /* end-point associated with f */
766 register ajint pi;
767 register ajint pj; /* end-point associated with p */
768 ajint cflag;
769 ajint rflag; /* for recomputation */
770 ajint *va; /* pointer to vv(A[i], B[j]) */
771 ajint limit; /* the bound on j */
772
773 /*
774 ** Reverse pass
775 ** rows
776 ** CC : the scores on the current row
777 ** RR and EE : the endpoints that lead to CC
778 ** DD : the deletion scores
779 ** SS and FF : the endpoints that lead to DD
780 **
781 ** columns
782 ** HH : the scores on the current columns
783 ** II and JJ : the endpoints that lead to HH
784 ** WW : the deletion scores
785 ** XX and YY : the endpoints that lead to WW
786 */
787
788 for(j = nn; j >= n1; j--)
789 {
790 CC[j] = 0;
791 EE[j] = j;
792 DD[j] = - (q);
793 FF[j] = j;
794 if(nseq == 2 || j > mm)
795 RR[j] = SS[j] = mm + 1;
796 else
797 RR[j] = SS[j] = j;
798 }
799
800 for(i = mm; i >= m1; i--)
801 {
802 c = p = 0;
803 f = - (q);
804 ci = fi = i;
805 pi = i + 1;
806 cj = fj = pj = nn + 1;
807 va = sub[(ajint)A[i]];
808
809 if(nseq == 2 || n1 > i)
810 limit = n1;
811 else
812 limit = i + 1;
813
814 for(j = nn; j >= limit; j--)
815 {
816 f = f - r;
817 c = c - qr;
818 MATCHERORDER(f, fi, fj, c, ci, cj)
819 c = CC[j] - qr;
820 ci = RR[j];
821 cj = EE[j];
822 d = DD[j] - r;
823 di = SS[j];
824 dj = FF[j];
825 MATCHERORDER(d, di, dj, c, ci, cj)
826 c = 0;
827 MATCHERDIAG(i, j, c, p+va[(ajint)B[j]]) /* diagonal */
828
829 if(c <= 0)
830 {
831 c = 0;
832 ci = i;
833 cj = j;
834 }
835 else
836 {
837 ci = pi;
838 cj = pj;
839 }
840
841 MATCHERORDER(c, ci, cj, d, di, dj)
842 MATCHERORDER(c, ci, cj, f, fi, fj)
843 p = CC[j];
844 CC[j] = c;
845 pi = RR[j];
846 pj = EE[j];
847 RR[j] = ci;
848 EE[j] = cj;
849 DD[j] = d;
850 SS[j] = di;
851 FF[j] = dj;
852 if(c > lmin)
853 flag = 1;
854 }
855
856 if(nseq == 2 || i < n1)
857 {
858 HH[i] = CC[n1];
859 II[i] = RR[n1];
860 JJ[i] = EE[n1];
861 WW[i] = DD[n1];
862 XX[i] = SS[n1];
863 YY[i] = FF[n1];
864 }
865 }
866
867 for(rl = m1, cl = n1; ;)
868 {
869 for(rflag = cflag = 1;( rflag && m1 > 1 ) || ( cflag && n1 > 1 );)
870 {
871 if(rflag && m1 > 1) /* Compute one row */
872 {
873 rflag = 0;
874 m1--;
875 c = p = 0;
876 f = - (q);
877 ci = fi = m1;
878 pi = m1 + 1;
879 cj = fj = pj = nn + 1;
880 va = sub[(ajint)A[m1]];
881
882 for(j = nn; j >= n1; j--)
883 {
884 f = f - r;
885 c = c - qr;
886 MATCHERORDER(f, fi, fj, c, ci, cj)
887 c = CC[j] - qr;
888 ci = RR[j];
889 cj = EE[j];
890 d = DD[j] - r;
891 di = SS[j];
892 dj = FF[j];
893 MATCHERORDER(d, di, dj, c, ci, cj)
894 c = 0;
895 MATCHERDIAG(m1, j, c, p+va[(ajint)B[j]]) /* diagonal */
896
897 if(c <= 0)
898 {
899 c = 0;
900 ci = m1;
901 cj = j;
902 }
903 else
904 {
905 ci = pi;
906 cj = pj;
907 }
908
909 MATCHERORDER(c, ci, cj, d, di, dj)
910 MATCHERORDER(c, ci, cj, f, fi, fj)
911 p = CC[j];
912 CC[j] = c;
913 pi = RR[j];
914 pj = EE[j];
915 RR[j] = ci;
916 EE[j] = cj;
917 DD[j] = d;
918 SS[j] = di;
919 FF[j] = dj;
920 if(c > lmin)
921 flag = 1;
922 if(! rflag && ( (ci > rl && cj > cl) || (di > rl && dj >
923 cl)
924 || (fi > rl && fj > cl )))
925 rflag = 1;
926 }
927
928 HH[m1] = CC[n1];
929 II[m1] = RR[n1];
930 JJ[m1] = EE[n1];
931 WW[m1] = DD[n1];
932 XX[m1] = SS[n1];
933 YY[m1] = FF[n1];
934 if(! cflag && ( (ci > rl && cj > cl) || (di > rl && dj > cl)
935 || (fi > rl && fj > cl)))
936 cflag = 1;
937 }
938
939 if(nseq == 1 && n1 == (m1 + 1) && ! rflag)
940 cflag = 0;
941
942 if(cflag && n1 > 1) /* Compute one column */
943 {
944 cflag = 0;
945 n1--;
946 c = 0;
947 f = - (q);
948 cj = fj = n1;
949 va = sub[(ajint)B[n1]];
950 if(nseq == 2 || mm < n1)
951 {
952 p = 0;
953 ci = fi = pi = mm + 1;
954 pj = n1 + 1;
955 limit = mm;
956 }
957 else
958 {
959 p = HH[n1];
960 pi = II[n1];
961 pj = JJ[n1];
962 ci = fi = n1;
963 limit = n1 - 1;
964 }
965
966 for(i = limit; i >= m1; i--)
967 {
968 f = f - r;
969 c = c - qr;
970 MATCHERORDER(f, fi, fj, c, ci, cj)
971 c = HH[i] - qr;
972 ci = II[i];
973 cj = JJ[i];
974 d = WW[i] - r;
975 di = XX[i];
976 dj = YY[i];
977 MATCHERORDER(d, di, dj, c, ci, cj)
978 c = 0;
979 MATCHERDIAG(i, n1, c, p+va[(ajint)A[i]])
980 if(c <= 0)
981 {
982 c = 0;
983 ci = i;
984 cj = n1;
985 }
986 else
987 {
988 ci = pi;
989 cj = pj;
990 }
991
992 MATCHERORDER(c, ci, cj, d, di, dj)
993 MATCHERORDER(c, ci, cj, f, fi, fj)
994 p = HH[i];
995 HH[i] = c;
996 pi = II[i];
997 pj = JJ[i];
998 II[i] = ci;
999 JJ[i] = cj;
1000 WW[i] = d;
1001 XX[i] = di;
1002 YY[i] = dj;
1003
1004 if(c > lmin)
1005 flag = 1;
1006
1007 if(! cflag && ( (ci > rl && cj > cl) || (di > rl && dj >
1008 cl)
1009 || (fi > rl && fj > cl)))
1010 cflag = 1;
1011 }
1012 CC[n1] = HH[m1];
1013 RR[n1] = II[m1];
1014 EE[n1] = JJ[m1];
1015 DD[n1] = WW[m1];
1016 SS[n1] = XX[m1];
1017 FF[n1] = YY[m1];
1018
1019 if(! rflag && ( (ci > rl && cj > cl) || (di > rl && dj > cl)
1020 || (fi > rl && fj > cl)))
1021 rflag = 1;
1022 }
1023 }
1024 if((m1 == 1 && n1 == 1) || matcher_NoCross(numnode))
1025 break;
1026 }
1027 m1--;
1028 n1--;
1029
1030 return 0;
1031 }
1032
1033
1034
1035
1036 /* @funcstatic matcher_SmallPass **********************************************
1037 **
1038 ** Undocumented
1039 **
1040 ** @param [r] A [const char*] Undocumented
1041 ** @param [r] B [const char*] Undocumented
1042 ** @param [r] count [ajint] Undocumented
1043 ** @param [r] nseq [ajint] Number of sequences
1044 ** @param [w] numnode [ajint*] Number of alignments in list
1045 ** @return [ajint] Undocumented
1046 ******************************************************************************/
1047
matcher_SmallPass(const char * A,const char * B,ajint count,ajint nseq,ajint * numnode)1048 static ajint matcher_SmallPass(const char *A,const char *B,
1049 ajint count,ajint nseq, ajint *numnode)
1050 {
1051 register ajint i;
1052 register ajint j; /* row and column indices */
1053 register ajint c; /* best score at current point */
1054 register ajint f; /* best score ending with insertion */
1055 register ajint d; /* best score ending with deletion */
1056 register ajint p; /* best score at (i-1, j-1) */
1057 register ajint ci;
1058 register ajint cj; /* end-point associated with c */
1059 register ajint di;
1060 register ajint dj; /* end-point associated with d */
1061 register ajint fi;
1062 register ajint fj; /* end-point associated with f */
1063 register ajint pi;
1064 register ajint pj; /* end-point associated with p */
1065 ajint *va; /* pointer to vv(A[i], B[j]) */
1066 ajint limit; /* lower bound on j */
1067
1068 for(j = n1 + 1; j <= nn; j++)
1069 {
1070 CC[j] = 0;
1071 RR[j] = m1;
1072 EE[j] = j;
1073 DD[j] = - (q);
1074 SS[j] = m1;
1075 FF[j] = j;
1076 }
1077
1078 for(i = m1 + 1; i <= mm; i++)
1079 {
1080 c = 0; /* Initialize column 0 */
1081 f = - (q);
1082 ci = fi = i;
1083 va = sub[(ajint)A[i]];
1084
1085 if(nseq == 2 || i <= n1)
1086 {
1087 p = 0;
1088 pi = i - 1;
1089 cj = fj = pj = n1;
1090 limit = n1 + 1;
1091 }
1092 else
1093 {
1094 p = CC[i];
1095 pi = RR[i];
1096 pj = EE[i];
1097 cj = fj = i;
1098 limit = i + 1;
1099 }
1100
1101 for(j = limit; j <= nn; j++)
1102 {
1103 f = f - r;
1104 c = c - qr;
1105 MATCHERORDER(f, fi, fj, c, ci, cj)
1106 c = CC[j] - qr;
1107 ci = RR[j];
1108 cj = EE[j];
1109 d = DD[j] - r;
1110 di = SS[j];
1111 dj = FF[j];
1112 MATCHERORDER(d, di, dj, c, ci, cj)
1113 c = 0;
1114 MATCHERDIAG(i, j, c, p+va[(ajint)B[j]]) /* diagonal */
1115
1116 if(c <= 0)
1117 {
1118 c = 0;
1119 ci = i;
1120 cj = j;
1121 }
1122 else
1123 {
1124 ci = pi;
1125 cj = pj;
1126 }
1127
1128 MATCHERORDER(c, ci, cj, d, di, dj)
1129 MATCHERORDER(c, ci, cj, f, fi, fj)
1130 p = CC[j];
1131 CC[j] = c;
1132 pi = RR[j];
1133 pj = EE[j];
1134 RR[j] = ci;
1135 EE[j] = cj;
1136 DD[j] = d;
1137 SS[j] = di;
1138 FF[j] = dj;
1139 if(c > lmin) /* add the score into list */
1140 lmin = matcher_Addnode(c, ci, cj, i, j, count, lmin, numnode);
1141 }
1142 }
1143
1144 return 0;
1145 }
1146
1147
1148
1149
1150 /* @funcstatic matcher_Addnode ************************************************
1151 **
1152 ** Add a node
1153 **
1154 ** @param [r] c [ajint] Undocumented
1155 ** @param [r] ci [ajint] Undocumented
1156 ** @param [r] cj [ajint] Undocumented
1157 ** @param [r] i [ajint] Undocumented
1158 ** @param [r] j [ajint] Undocumented
1159 ** @param [r] K [ajint] Undocumented
1160 ** @param [r] cost [ajint] Undocumented
1161 ** @param [w] numnode [ajint*] NUmber of nodes (updated)
1162 ** @return [ajint] Undocumented
1163 ******************************************************************************/
1164
matcher_Addnode(ajint c,ajint ci,ajint cj,ajint i,ajint j,ajint K,ajint cost,ajint * numnode)1165 static ajint matcher_Addnode(ajint c, ajint ci, ajint cj, ajint i, ajint j,
1166 ajint K, ajint cost, ajint *numnode)
1167 {
1168 ajint found; /* 1 if the node is in LIST */
1169 register ajint d;
1170
1171 found = 0;
1172
1173 if(most != 0 && most->STARI == ci && most->STARJ == cj)
1174 found = 1;
1175 else
1176 for(d = 0; d < *numnode; d++)
1177 {
1178 most = LIST[d];
1179 if(most->STARI == ci && most->STARJ == cj)
1180 {
1181 found = 1;
1182 break;
1183 }
1184 }
1185
1186 if(found)
1187 {
1188 if(most->SCORE < c)
1189 {
1190 most->SCORE = c;
1191 most->ENDI = i;
1192 most->ENDJ = j;
1193 }
1194
1195 if(most->TOP > i)
1196 most->TOP = i;
1197
1198 if(most->BOT < i)
1199 most->BOT = i;
1200
1201 if(most->LEFT > j)
1202 most->LEFT = j;
1203
1204 if(most->RIGHT < j)
1205 most->RIGHT = j;
1206 }
1207 else
1208 {
1209 if(*numnode == K) /* list full */
1210 most = low;
1211 else
1212 most = LIST[(*numnode)++];
1213 most->SCORE = c;
1214 most->STARI = ci;
1215 most->STARJ = cj;
1216 most->ENDI = i;
1217 most->ENDJ = j;
1218 most->TOP = most->BOT = i;
1219 most->LEFT = most->RIGHT = j;
1220 }
1221
1222 if(*numnode == K)
1223 {
1224 if(low == most || ! low)
1225 for(low = LIST[0], d = 1; d < *numnode; d++)
1226 if(LIST[d]->SCORE < low->SCORE)
1227 low = LIST[d];
1228
1229 return (low->SCORE);
1230 }
1231
1232
1233 return cost;
1234 }
1235
1236
1237
1238
1239 /* @funcstatic matcher_Findmax ************************************************
1240 **
1241 ** Find maximum node
1242 **
1243 ** @param [u] numnode [ajint*] Maximum node number
1244 ** @return [vertexptr] Maximum node
1245 ******************************************************************************/
1246
matcher_Findmax(ajint * numnode)1247 static vertexptr matcher_Findmax(ajint *numnode)
1248 {
1249 vertexptr cur;
1250 register ajint i;
1251 register ajint j;
1252
1253 for(j = 0, i = 1; i < *numnode; i++)
1254 if(LIST[i]->SCORE > LIST[j]->SCORE)
1255 j = i;
1256 cur = LIST[j];
1257
1258 if(j != --(*numnode))
1259 {
1260 LIST[j] = LIST[*numnode];
1261 LIST[*numnode] = cur;
1262 }
1263 most = LIST[0];
1264 if(low == cur)
1265 low = LIST[0];
1266
1267 return (cur);
1268 }
1269
1270
1271
1272
1273 /* @funcstatic matcher_Diff ***************************************************
1274 **
1275 ** Undocumented
1276 **
1277 ** @param [r] A [const char*] Undocumented
1278 ** @param [r] B [const char*] Undocumented
1279 ** @param [r] M [ajint] Undocumented
1280 ** @param [r] N [ajint] Undocumented
1281 ** @param [r] tb [ajint] Undocumented
1282 ** @param [r] te [ajint] Undocumented
1283 ** @return [ajint] Undocumented
1284 ******************************************************************************/
1285
1286
matcher_Diff(const char * A,const char * B,ajint M,ajint N,ajint tb,ajint te)1287 static ajint matcher_Diff(const char *A,const char *B,
1288 ajint M,ajint N,ajint tb,ajint te)
1289 {
1290 ajint midi;
1291 ajint midj;
1292 ajint type; /* Midpoint, type, and cost */
1293 ajint midc;
1294 ajint zero = 0; /* ajint type zero */
1295
1296 register ajint i;
1297 register ajint j;
1298 register ajint c;
1299 register ajint e;
1300 register ajint d;
1301 register ajint s;
1302 ajint t;
1303 ajint *va;
1304
1305 /* Boundary cases: M <= 1 or N == 0 */
1306
1307 if(N <= 0)
1308 {
1309 if(M > 0)
1310 MATCHERDEL(M)
1311 return - matchergap(M);
1312 }
1313
1314 if(M <= 1)
1315 {
1316 if(M <= 0)
1317 {
1318 MATCHERINS(N);
1319 return - matchergap(N);
1320 }
1321 if(tb > te)
1322 tb = te;
1323 midc = - (tb + r + matchergap(N));
1324 midj = 0;
1325 va = sub[(ajint)A[1]];
1326
1327 for(j = 1; j <= N; j++)
1328 {
1329 for(tt = 1, z = row[I+1]; z != PAIRNULL; z = z->NEXT)
1330 if(z->COL == j+J)
1331 {
1332 tt = 0;
1333 break;
1334 }
1335
1336 if(tt)
1337 {
1338 c = va[(ajint)B[j]] - ( matchergap(j-1) + matchergap(N-j));
1339 if(c > midc)
1340 {
1341 midc = c;
1342 midj = j;
1343 }
1344 }
1345 }
1346
1347 if(midj == 0)
1348 {
1349 MATCHERINS(N)
1350 MATCHERDEL(1)
1351 }
1352 else
1353 {
1354 if(midj > 1) MATCHERINS(midj-1)
1355 MATCHERREP
1356
1357 /* mark(A[I],B[J]) as used: put J into list row[I] */
1358 I++; J++;
1359 AJNEW0(z);
1360 z->COL = J;
1361 z->NEXT = row[I];
1362 row[I] = z;
1363
1364 if(midj < N)
1365 MATCHERINS(N-midj)
1366 }
1367 return midc;
1368 }
1369
1370 /* Divide: Find optimum midpoint (midi,midj) of cost midc */
1371
1372 midi = M/2; /* Forward phase: */
1373 CC[0] = 0; /* Compute C(M/2,k) & D(M/2,k) for all k */
1374
1375 t = -q;
1376 for(j = 1; j <= N; j++)
1377 {
1378 CC[j] = t = t-r;
1379 DD[j] = t-q;
1380 }
1381 t = -tb;
1382
1383 for(i = 1; i <= midi; i++)
1384 {
1385 s = CC[0];
1386 CC[0] = c = t = t-r;
1387 e = t-q;
1388 va = sub[(ajint)A[i]];
1389 for(j = 1; j <= N; j++)
1390 {
1391 if((c = c - qr) > (e = e - r))
1392 e = c;
1393
1394 if((c = CC[j] - qr) > (d = DD[j] - r))
1395 d = c;
1396 MATCHERDIAG(i+I, j+J, c, s+va[(ajint)B[j]])
1397
1398 if(c < d)
1399 c = d;
1400
1401 if(c < e)
1402 c = e;
1403
1404 s = CC[j];
1405 CC[j] = c;
1406 DD[j] = d;
1407 }
1408 }
1409 DD[0] = CC[0];
1410
1411 RR[N] = 0; /* Reverse phase: */
1412 t = -q; /* Compute R(M/2,k) & S(M/2,k) for all k */
1413
1414 for(j = N-1; j >= 0; j--)
1415 {
1416 RR[j] = t = t-r;
1417 SS[j] = t-q;
1418 }
1419 t = -te;
1420
1421 for(i = M-1; i >= midi; i--)
1422 {
1423 s = RR[N];
1424 RR[N] = c = t = t-r;
1425 e = t-q;
1426 va = sub[(ajint)A[i+1]];
1427
1428 for(j = N-1; j >= 0; j--)
1429 {
1430 if((c = c - qr) > (e = e - r))
1431 e = c;
1432
1433 if((c = RR[j] - qr) > (d = SS[j] - r))
1434 d = c;
1435 MATCHERDIAG(i+1+I, j+1+J, c, s+va[(ajint)B[j+1]])
1436
1437 if(c < d)
1438 c = d;
1439
1440 if(c < e)
1441 c = e;
1442
1443 s = RR[j];
1444 RR[j] = c;
1445 SS[j] = d;
1446 }
1447 }
1448 SS[N] = RR[N];
1449
1450 midc = CC[0]+RR[0]; /* Find optimal midpoint */
1451 midj = 0;
1452 type = 1;
1453
1454 for(j = 0; j <= N; j++)
1455 if((c = CC[j] + RR[j]) >= midc)
1456 if(c > midc || (CC[j] != DD[j] && RR[j] == SS[j]))
1457 {
1458 midc = c;
1459 midj = j;
1460 }
1461
1462 for(j = N; j >= 0; j--)
1463 if((c = DD[j] + SS[j] + q) > midc)
1464 {
1465 midc = c;
1466 midj = j;
1467 type = 2;
1468 }
1469
1470
1471 /* Conquer: recursively around midpoint */
1472
1473 if(type == 1)
1474 {
1475 matcher_Diff(A,B,midi,midj,tb,q);
1476 matcher_Diff(A+midi,B+midj,M-midi,N-midj,q,te);
1477 }
1478 else
1479 {
1480 matcher_Diff(A,B,midi-1,midj,tb,zero);
1481 MATCHERDEL(2);
1482 matcher_Diff(A+midi+1,B+midj,M-midi-1,N-midj,zero,te);
1483 }
1484
1485 return midc;
1486 }
1487
1488
1489
1490
1491 /* @funcstatic matcher_Calcons ************************************************
1492 **
1493 ** Calculate a consensus sequence
1494 **
1495 ** @param [w] seqc0 [char*] Sequence A alignment output
1496 ** @param [w] seqc1 [char*] Sequence B alignment output
1497 ** @param [r] res [const ajint*] Undocumented
1498 ** @param [r] min0 [ajint] Undocumented
1499 ** @param [r] min1 [ajint] Undocumented
1500 ** @param [r] max0 [ajint] Undocumented
1501 ** @param [r] max1 [ajint] Undocumented
1502 ** @param [w] nc [ajint*] Number of conserved positions
1503 ** @param [w] nident [ajint*] Number of identical positions
1504 ** @return [ajint] Undocumented
1505 ******************************************************************************/
1506
matcher_Calcons(char * seqc0,char * seqc1,const ajint * res,ajint min0,ajint min1,ajint max0,ajint max1,ajint * nc,ajint * nident)1507 static ajint matcher_Calcons(char *seqc0,char *seqc1,
1508 const ajint *res,
1509 ajint min0, ajint min1,
1510 ajint max0, ajint max1,
1511 ajint *nc,ajint *nident)
1512 {
1513 ajint i0;
1514 ajint i1;
1515 ajint op;
1516 ajint nid;
1517 ajint lenc;
1518 ajint nd;
1519 char *sp0;
1520 char *sp1;
1521 const ajint *rp;
1522
1523 const char *sq1;
1524 const char *sq2;
1525
1526 /* now get the middle */
1527
1528 sp0 = seqc0 /*+mins*/;
1529 sp1 = seqc1 /*+mins*/;
1530 rp = res;
1531 lenc = nid = op = 0;
1532 i0 = min0;
1533 i1 = min1;
1534
1535 sq1 = ajStrGetPtr(ajSeqGetSeqS(seq));
1536 sq2 = ajStrGetPtr(ajSeqGetSeqS(seq2));
1537
1538 while(i0 < max0 || i1 < max1)
1539 {
1540 if(op == 0 && *rp == 0)
1541 {
1542 op = *rp++;
1543 *sp0 = sq1[i0++];
1544 *sp1 = sq2[i1++];
1545 lenc++;
1546 if(*sp0 == *sp1)
1547 nid++;
1548
1549 sp0++; sp1++;
1550 }
1551 else
1552 {
1553 if(op==0)
1554 op = *rp++;
1555
1556 if(op>0)
1557 {
1558 *sp0++ = '-';
1559 *sp1++ = sq2[i1++];
1560 op--;
1561 lenc++;
1562 }
1563 else
1564 {
1565 *sp0++ = sq1[i0++];
1566 *sp1++ = '-';
1567 op++;
1568 lenc++;
1569 }
1570 }
1571 }
1572 *sp0 = '\0';
1573 *sp1 = '\0';
1574 *nident = nid;
1575 *nc = lenc;
1576 /* get the right end */
1577 nd = 0;
1578
1579 return lenc+nd;
1580 }
1581