1 #include "mltaln.h"
2
3 #if 0
4 static FILE *fftfp;
5 #endif
6 static TLS int n20or4or2;
7
8 #define KEIKA 0
9 #define RND 0
10 #define DEBUG 0
11
12 #if RND // by D.Mathog
generateRndSeq(char * seq,int len)13 static void generateRndSeq( char *seq, int len )
14 {
15 while( len-- )
16 #if 1
17 *seq++ = (int)( rnd() * n20or4or2 );
18 #else
19 *seq++ = (int)1;
20 #endif
21 }
22 #endif
23
vec_init(Fukusosuu * result,int nlen)24 static void vec_init( Fukusosuu *result, int nlen )
25 {
26 while( nlen-- )
27 {
28 result->R = result->I = 0.0;
29 result++;
30 }
31 }
32
33 #if 0 // by D.Mathog
34 static void vec_init2( Fukusosuu **result, char *seq, double eff, int st, int ed )
35 {
36 int i;
37 for( i=st; i<ed; i++ )
38 result[(int)*seq++][i].R += eff;
39 }
40 #endif
41
seq_vec_2(Fukusosuu * result,double * score,double incr,char * seq)42 static void seq_vec_2( Fukusosuu *result, double *score, double incr, char *seq )
43 {
44 static TLS int n;
45 for( ; *seq; result++ )
46 {
47 n = amino_n[(int)*seq++];
48 if( n < 20 && n >= 0 ) result->R += incr * score[n];
49 #if 0
50 fprintf( stderr, "n=%d, score=%f, inc=%f R=%f\n",n, score[n], incr * score[n], result->R );
51 #endif
52 }
53 }
54
seq_vec_3(Fukusosuu ** result,double incr,char * seq)55 static void seq_vec_3( Fukusosuu **result, double incr, char *seq )
56 {
57 int i;
58 int n;
59 for( i=0; *seq; i++ )
60 {
61 n = amino_n[(int)*seq++];
62 if( n < n20or4or2 && n >= 0 ) result[n][i].R += incr;
63 }
64 }
65
seq_vec_5(Fukusosuu * result,double * score1,double * score2,double incr,char * seq)66 static void seq_vec_5( Fukusosuu *result, double *score1, double *score2, double incr, char *seq )
67 {
68 int n;
69 for( ; *seq; result++ )
70 {
71 n = amino_n[(int)*seq++];
72 if( n > 20 ) continue;
73 result->R += incr * score1[n];
74 result->I += incr * score2[n];
75 #if 0
76 fprintf( stderr, "n=%d, score=%f, inc=%f R=%f\n",n, score[n], incr * score[n], result->R );
77 #endif
78 }
79 }
80
81
seq_vec_4(Fukusosuu * result,double incr,char * seq)82 static void seq_vec_4( Fukusosuu *result, double incr, char *seq )
83 {
84 char s;
85 for( ; *seq; result++ )
86 {
87 s = *seq++;
88 if( s == 'a' )
89 result->R += incr;
90 else if( s == 't' )
91 result->R -= incr;
92 else if( s == 'g' )
93 result->I += incr;
94 else if( s == 'c' )
95 result->I -= incr;
96 }
97 }
98
99 #if 0 // by D.Mathog
100 static void seq_vec( Fukusosuu *result, char query, double incr, char *seq )
101 {
102 #if 0
103 int bk = nlen;
104 #endif
105 while( *seq )
106 {
107 if( *seq++ == query ) result->R += incr;
108 result++;
109 #if 0
110 fprintf( stderr, "i = %d result->R = %f\n", bk-nlen, (result-1)->R );
111 #endif
112 }
113 }
114
115 static int checkRepeat( int num, int *cutpos )
116 {
117 int tmp, buf;
118
119 buf = *cutpos;
120 while( num-- )
121 {
122 if( ( tmp = *cutpos++ ) < buf ) return( 1 );
123 buf = tmp;
124 }
125 return( 0 );
126 }
127
128 static int segcmp( void *ptr1, void *ptr2 )
129 {
130 int diff;
131 Segment **seg1 = (Segment **)ptr1;
132 Segment **seg2 = (Segment **)ptr2;
133 #if 0
134 return( (*seg1)->center - (*seg2)->center );
135 #else
136 diff = (*seg1)->center - (*seg2)->center;
137 if( diff ) return( diff );
138
139 diff = (*seg1)->start - (*seg2)->start;
140 if( diff ) return( diff );
141
142 diff = (*seg1)->end - (*seg2)->end;
143 if( diff ) return( diff );
144
145 fprintf( stderr, "USE STABLE SORT !!\n" );
146 exit( 1 );
147 return( 0 );
148 #endif
149 }
150 #endif
151
152
mymergesort(int first,int last,Segment ** seg)153 static void mymergesort( int first, int last, Segment **seg )
154 {
155 int middle;
156 static TLS int i, j, k, p;
157 static TLS int allo = 0;
158 static TLS Segment **work = NULL;
159
160 if( seg == NULL )
161 {
162 if( work ) free( work );
163 work = NULL;
164 allo = 0;
165 return;
166 }
167
168 if( last > allo )
169 {
170 allo = last;
171 if( work ) free( work );
172 work = (Segment **)calloc( allo / 2 + 1, sizeof( Segment *) );
173 }
174
175 if( first < last )
176 {
177 middle = ( first + last ) / 2;
178 mymergesort( first, middle, seg );
179 mymergesort( middle+1, last, seg );
180 p = 0;
181 for( i=first; i<=middle; i++ ) work[p++] = seg[i];
182 i = middle + 1; j = 0; k = first;
183 while( i <= last && j < p )
184 {
185 if( work[j]->center <= seg[i]->center )
186 seg[k++] = work[j++];
187 else
188 seg[k++] = seg[i++];
189 }
190 while( j < p ) seg[k++] = work[j++];
191 }
192 }
193
194
Fgetlag(double ** n_dynamicmtx,char ** seq1,char ** seq2,double * eff1,double * eff2,int clus1,int clus2,int alloclen)195 double Fgetlag(
196 double **n_dynamicmtx,
197 char **seq1, char **seq2,
198 double *eff1, double *eff2,
199 int clus1, int clus2,
200 int alloclen )
201 {
202 int i, j, k, l, m;
203 int nlen, nlen2, nlen4;
204 static TLS int crossscoresize = 0;
205 static TLS char **tmpseq1 = NULL;
206 static TLS char **tmpseq2 = NULL;
207 static TLS char **tmpptr1 = NULL;
208 static TLS char **tmpptr2 = NULL;
209 static TLS char **tmpres1 = NULL;
210 static TLS char **tmpres2 = NULL;
211 static TLS char **result1 = NULL;
212 static TLS char **result2 = NULL;
213 #if RND
214 static TLS char **rndseq1 = NULL;
215 static TLS char **rndseq2 = NULL;
216 #endif
217 static TLS Fukusosuu **seqVector1 = NULL;
218 static TLS Fukusosuu **seqVector2 = NULL;
219 static TLS Fukusosuu **naiseki = NULL;
220 static TLS Fukusosuu *naisekiNoWa = NULL;
221 static TLS double *soukan = NULL;
222 static TLS double **crossscore = NULL;
223 int nlentmp;
224 static TLS int *kouho = NULL;
225 static TLS Segment *segment = NULL;
226 static TLS Segment *segment1 = NULL;
227 static TLS Segment *segment2 = NULL;
228 static TLS Segment **sortedseg1 = NULL;
229 static TLS Segment **sortedseg2 = NULL;
230 static TLS int *cut1 = NULL;
231 static TLS int *cut2 = NULL;
232 static TLS int localalloclen = 0;
233 int lag;
234 int tmpint;
235 int count, count0;
236 int len1, len2;
237 int totallen;
238 double dumdb = 0.0;
239 int headgp, tailgp;
240
241 len1 = strlen( seq1[0] );
242 len2 = strlen( seq2[0] );
243 nlentmp = MAX( len1, len2 );
244
245 nlen = 1;
246 while( nlentmp >= nlen ) nlen <<= 1;
247 #if 0
248 fprintf( stderr, "### nlen = %d\n", nlen );
249 #endif
250
251 nlen2 = nlen/2; nlen4 = nlen2 / 2;
252
253 #if DEBUG
254 fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
255 fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
256 #endif
257
258 if( !localalloclen )
259 {
260 kouho = AllocateIntVec( NKOUHO );
261 cut1 = AllocateIntVec( MAXSEG );
262 cut2 = AllocateIntVec( MAXSEG );
263 tmpptr1 = AllocateCharMtx( njob, 0 );
264 tmpptr2 = AllocateCharMtx( njob, 0 );
265 result1 = AllocateCharMtx( njob, alloclen );
266 result2 = AllocateCharMtx( njob, alloclen );
267 tmpres1 = AllocateCharMtx( njob, alloclen );
268 tmpres2 = AllocateCharMtx( njob, alloclen );
269 // crossscore = AllocateDoubleMtx( MAXSEG, MAXSEG );
270 segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
271 segment1 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
272 segment2 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
273 sortedseg1 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
274 sortedseg2 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
275 if( !( segment && segment1 && segment2 && sortedseg1 && sortedseg2 ) )
276 ErrorExit( "Allocation error\n" );
277
278 if ( scoremtx == -1 ) n20or4or2 = 4;
279 else if( fftscore == 1 ) n20or4or2 = 2;
280 else n20or4or2 = 20;
281 }
282 if( localalloclen < nlen )
283 {
284 if( localalloclen )
285 {
286 #if 1
287 FreeFukusosuuMtx ( seqVector1 );
288 FreeFukusosuuMtx ( seqVector2 );
289 FreeFukusosuuVec( naisekiNoWa );
290 FreeFukusosuuMtx( naiseki );
291 FreeDoubleVec( soukan );
292 FreeCharMtx( tmpseq1 );
293 FreeCharMtx( tmpseq2 );
294 #endif
295 #if RND
296 FreeCharMtx( rndseq1 );
297 FreeCharMtx( rndseq2 );
298 #endif
299 }
300
301
302 tmpseq1 = AllocateCharMtx( njob, nlen );
303 tmpseq2 = AllocateCharMtx( njob, nlen );
304 naisekiNoWa = AllocateFukusosuuVec( nlen );
305 naiseki = AllocateFukusosuuMtx( n20or4or2, nlen );
306 seqVector1 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
307 seqVector2 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
308 soukan = AllocateDoubleVec( nlen+1 );
309
310 #if RND
311 rndseq1 = AllocateCharMtx( njob, nlen );
312 rndseq2 = AllocateCharMtx( njob, nlen );
313 for( i=0; i<njob; i++ )
314 {
315 generateRndSeq( rndseq1[i], nlen );
316 generateRndSeq( rndseq2[i], nlen );
317 }
318 #endif
319 localalloclen = nlen;
320 }
321
322 for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
323 for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );
324
325 #if 0
326 fftfp = fopen( "input_of_Falign", "w" );
327 fprintf( fftfp, "nlen = %d\n", nlen );
328 fprintf( fftfp, "seq1: ( %d sequences ) \n", clus1 );
329 for( i=0; i<clus1; i++ )
330 fprintf( fftfp, "%s\n", seq1[i] );
331 fprintf( fftfp, "seq2: ( %d sequences ) \n", clus2 );
332 for( i=0; i<clus2; i++ )
333 fprintf( fftfp, "%s\n", seq2[i] );
334 fclose( fftfp );
335 system( "less input_of_Falign < /dev/tty > /dev/tty" );
336 #endif
337
338 if( fftkeika ) fprintf( stderr, " FFT ... " );
339
340 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
341 if( fftscore && scoremtx != -1 )
342 {
343 for( i=0; i<clus1; i++ )
344 {
345 seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
346 seq_vec_2( seqVector1[1], volume, eff1[i], tmpseq1[i] );
347 }
348 }
349 else
350 {
351 #if 0
352 for( i=0; i<clus1; i++ ) for( j=0; j<n20or4or2; j++ )
353 seq_vec( seqVector1[j], amino[j], eff1[i], tmpseq1[i] );
354 #else
355 for( i=0; i<clus1; i++ )
356 seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
357 #endif
358 }
359 #if RND
360 for( i=0; i<clus1; i++ )
361 {
362 vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
363 }
364 #endif
365 #if 0
366 fftfp = fopen( "seqVec", "w" );
367 fprintf( fftfp, "before transform\n" );
368 for( k=0; k<n20or4or2; k++ )
369 {
370 fprintf( fftfp, "nlen=%d\n", nlen );
371 fprintf( fftfp, "%c\n", amino[k] );
372 for( l=0; l<nlen; l++ )
373 fprintf( fftfp, "%f %f\n", seqVector1[k][l].R, seqVector1[k][l].I );
374 }
375 fclose( fftfp );
376 system( "less seqVec < /dev/tty > /dev/tty" );
377 #endif
378
379 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
380 if( fftscore && scoremtx != -1 )
381 {
382 for( i=0; i<clus2; i++ )
383 {
384 seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
385 seq_vec_2( seqVector2[1], volume, eff2[i], tmpseq2[i] );
386 }
387 }
388 else
389 {
390 #if 0
391 for( i=0; i<clus2; i++ ) for( j=0; j<n20or4or2; j++ )
392 seq_vec( seqVector2[j], amino[j], eff2[i], tmpseq2[i] );
393 #else
394 for( i=0; i<clus2; i++ )
395 seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
396 #endif
397 }
398 #if RND
399 for( i=0; i<clus2; i++ )
400 {
401 vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
402 }
403 #endif
404
405 #if 0
406 fftfp = fopen( "seqVec2", "w" );
407 fprintf( fftfp, "before fft\n" );
408 for( k=0; k<n20or4or2; k++ )
409 {
410 fprintf( fftfp, "%c\n", amino[k] );
411 for( l=0; l<nlen; l++ )
412 fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
413 }
414 fclose( fftfp );
415 system( "less seqVec2 < /dev/tty > /dev/tty" );
416 #endif
417
418 for( j=0; j<n20or4or2; j++ )
419 {
420 fft( nlen, seqVector2[j], 0 );
421 fft( nlen, seqVector1[j], 0 );
422 }
423 #if 0
424 fftfp = fopen( "seqVec2", "w" );
425 fprintf( fftfp, "#after fft\n" );
426 for( k=0; k<n20or4or2; k++ )
427 {
428 fprintf( fftfp, "#%c\n", amino[k] );
429 for( l=0; l<nlen; l++ )
430 fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
431 }
432 fclose( fftfp );
433 system( "less seqVec2 < /dev/tty > /dev/tty" );
434 #endif
435
436 for( k=0; k<n20or4or2; k++ )
437 {
438 for( l=0; l<nlen; l++ )
439 calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
440 }
441 for( l=0; l<nlen; l++ )
442 {
443 naisekiNoWa[l].R = 0.0;
444 naisekiNoWa[l].I = 0.0;
445 for( k=0; k<n20or4or2; k++ )
446 {
447 naisekiNoWa[l].R += naiseki[k][l].R;
448 naisekiNoWa[l].I += naiseki[k][l].I;
449 }
450 }
451
452 #if 0
453 fftfp = fopen( "naisekiNoWa", "w" );
454 fprintf( fftfp, "#Before fft\n" );
455 for( l=0; l<nlen; l++ )
456 fprintf( fftfp, "%d %f %f\n", l, naisekiNoWa[l].R, naisekiNoWa[l].I );
457 fclose( fftfp );
458 system( "less naisekiNoWa < /dev/tty > /dev/tty " );
459 #endif
460
461 fft( -nlen, naisekiNoWa, 0 );
462
463 for( m=0; m<=nlen2; m++ )
464 soukan[m] = naisekiNoWa[nlen2-m].R;
465 for( m=nlen2+1; m<nlen; m++ )
466 soukan[m] = naisekiNoWa[nlen+nlen2-m].R;
467
468 #if 0
469 fftfp = fopen( "naisekiNoWa", "w" );
470 fprintf( fftfp, "#After fft\n" );
471 for( l=0; l<nlen; l++ )
472 fprintf( fftfp, "%d %f\n", l, naisekiNoWa[l].R );
473 fclose( fftfp );
474 fftfp = fopen( "list.plot", "w" );
475 fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
476 fclose( fftfp );
477 system( "/usr/bin/gnuplot list.plot &" );
478 #endif
479 #if 0
480 fprintf( stderr, "frt write start\n" );
481 fftfp = fopen( "frt", "w" );
482 for( l=0; l<nlen; l++ )
483 fprintf( fftfp, "%d %f\n", l-nlen2, soukan[l] );
484 fclose( fftfp );
485 system( "less frt < /dev/tty > /dev/tty" );
486 #if 0
487 fftfp = fopen( "list.plot", "w" );
488 fprintf( fftfp, "plot 'frt'\n pause +1" );
489 fclose( fftfp );
490 system( "/usr/bin/gnuplot list.plot" );
491 #endif
492 #endif
493
494
495 getKouho( kouho, NKOUHO, soukan, nlen );
496
497 #if 0
498 for( i=0; i<NKOUHO; i++ )
499 {
500 fprintf( stdout, "kouho[%d] = %d\n", i, kouho[i] );
501 }
502 #endif
503
504 #if KEIKA
505 fprintf( stderr, "Searching anchors ... " );
506 #endif
507 count = 0;
508
509
510
511 #define CAND 0
512 #if CAND
513 fftfp = fopen( "cand", "w" );
514 fclose( fftfp );
515 #endif
516
517 for( k=0; k<NKOUHO; k++ )
518 {
519
520 lag = kouho[k];
521 zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
522 #if CAND
523 fftfp = fopen( "cand", "a" );
524 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
525 fprintf( fftfp, "%s\n", tmpptr1[0] );
526 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
527 fprintf( fftfp, "%s\n", tmpptr2[0] );
528 fprintf( fftfp, ">\n", k+1, lag );
529 fclose( fftfp );
530 #endif
531 tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );
532
533 if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );
534
535
536 if( tmpint == 0 ) break; // 060430 iinoka ?
537 while( tmpint-- > 0 )
538 {
539 if( lag > 0 )
540 {
541 segment1[count].start = segment[count].start ;
542 segment1[count].end = segment[count].end ;
543 segment1[count].center = segment[count].center;
544 segment1[count].score = segment[count].score;
545
546 segment2[count].start = segment[count].start + lag;
547 segment2[count].end = segment[count].end + lag;
548 segment2[count].center = segment[count].center + lag;
549 segment2[count].score = segment[count].score ;
550 }
551 else
552 {
553 segment1[count].start = segment[count].start - lag;
554 segment1[count].end = segment[count].end - lag;
555 segment1[count].center = segment[count].center - lag;
556 segment1[count].score = segment[count].score ;
557
558 segment2[count].start = segment[count].start ;
559 segment2[count].end = segment[count].end ;
560 segment2[count].center = segment[count].center;
561 segment2[count].score = segment[count].score ;
562 }
563 #if 0
564 fprintf( stderr, "Goukaku=%dko\n", tmpint );
565 fprintf( stderr, "in 1 %d\n", segment1[count].center );
566 fprintf( stderr, "in 2 %d\n", segment2[count].center );
567 #endif
568 segment1[count].pair = &segment2[count];
569 segment2[count].pair = &segment1[count];
570 count++;
571 #if 0
572 fprintf( stderr, "count=%d\n", count );
573 #endif
574 }
575 }
576
577 #if 1
578 fprintf( stderr, "done. (%d anchors)\r", count );
579 #endif
580 if( !count && fftNoAnchStop )
581 ErrorExit( "Cannot detect anchor!" );
582 #if 0
583 fprintf( stdout, "RESULT before sort:\n" );
584 for( l=0; l<count+1; l++ )
585 {
586 fprintf( stdout, "cut[%d]=%d, ", l, segment1[l].center );
587 fprintf( stdout, "%d score = %f\n", segment2[l].center, segment1[l].score );
588 }
589 exit( 1 );
590 #endif
591
592 #if KEIKA
593 fprintf( stderr, "Aligning anchors ... " );
594 #endif
595 for( i=0; i<count; i++ )
596 {
597 sortedseg1[i] = &segment1[i];
598 sortedseg2[i] = &segment2[i];
599 }
600
601 {
602 mymergesort( 0, count-1, sortedseg1 );
603 mymergesort( 0, count-1, sortedseg2 );
604 for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
605 for( i=0; i<count; i++ ) sortedseg2[i]->number = i;
606
607 if( crossscoresize < count+2 )
608 {
609 crossscoresize = count+2;
610 fprintf( stderr, "####################################################################################################################################allocating crossscore, size = %d\n", crossscoresize );
611 if( crossscore ) FreeDoubleMtx( crossscore );
612 crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
613 }
614
615 for( i=0; i<count+2; i++ ) for( j=0; j<count+2; j++ )
616 crossscore[i][j] = 0.0;
617 for( i=0; i<count; i++ )
618 {
619 crossscore[segment1[i].number+1][segment1[i].pair->number+1] = segment1[i].score;
620 cut1[i+1] = sortedseg1[i]->center;
621 cut2[i+1] = sortedseg2[i]->center;
622 }
623
624 #if DEBUG
625 fprintf( stderr, "AFTER SORT\n" );
626 for( i=0; i<count; i++ ) fprintf( stderr, "%d, %d\n", segment1[i].start, segment2[i].start );
627 #endif
628
629 crossscore[0][0] = 10000000.0;
630 cut1[0] = 0;
631 cut2[0] = 0;
632 crossscore[count+1][count+1] = 10000000.0;
633 cut1[count+1] = len1;
634 cut2[count+1] = len2;
635 count += 2;
636 count0 = count;
637
638 blockAlign2( cut1, cut2, sortedseg1, sortedseg2, crossscore, &count );
639 }
640 if( fftkeika )
641 {
642 if( count0 > count )
643 {
644 fprintf( stderr, "REPEAT!? \n" );
645 if( fftRepeatStop ) exit( 1 );
646 }
647 #if KEIKA
648 else
649 fprintf( stderr, "done\n" );
650 fprintf( stderr, "done. (%d anchors)\n", count );
651 #endif
652 }
653
654 #if 0
655 fftfp = fopen( "fft", "a" );
656 fprintf( fftfp, "RESULT after sort:\n" );
657 for( l=0; l<count; l++ )
658 {
659 fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
660 fprintf( fftfp, "%d\n", segment2[l].center );
661 }
662 fclose( fftfp );
663 #endif
664
665 #if 0
666 fftfp = fopen( "fft", "a" );
667 fprintf( fftfp, "RESULT after sort:\n" );
668 for( l=0; l<count; l++ )
669 {
670 fprintf( fftfp, "cut : %d %d\n", cut1[l], cut2[l] );
671 }
672 fclose( fftfp );
673 #endif
674
675 #if KEIKA
676 fprintf( trap_g, "Devided to %d segments\n", count-1 );
677 fprintf( trap_g, "%d %d forg\n", MIN( clus1, clus2 ), count-1 );
678 #endif
679
680 totallen = 0;
681 for( j=0; j<clus1; j++ ) result1[j][0] = 0;
682 for( j=0; j<clus2; j++ ) result2[j][0] = 0;
683 for( i=0; i<count-1; i++ )
684 {
685 if( i == 0 ) headgp = outgap; else headgp = 1;
686 if( i == count-2 ) tailgp = outgap; else tailgp = 1;
687
688 #if DEBUG
689 fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
690 #else
691 #if KEIKA
692 fprintf( stderr, "DP %03d / %03d\r", i+1, count-1 );
693 #endif
694 #endif
695 for( j=0; j<clus1; j++ )
696 {
697 strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
698 tmpres1[j][cut1[i+1]-cut1[i]] = 0;
699 }
700 for( j=0; j<clus2; j++ )
701 {
702 strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
703 tmpres2[j][cut2[i+1]-cut2[i]] = 0;
704 }
705 switch( alg )
706 {
707 case( 'a' ):
708 Aalign( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen );
709 break;
710 case( 'M' ):
711 MSalignmm( n_dynamicmtx, tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, headgp, tailgp );
712 break;
713 case( 'A' ):
714 if( clus1 == 1 && clus2 == 1 )
715 G__align11( n_dynamicmtx, tmpres1, tmpres2, alloclen, headgp, tailgp );
716 else
717 A__align( n_dynamicmtx, tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumdb, NULL, NULL, NULL, NULL, NULL, 0, NULL, headgp, tailgp );
718 break;
719 default:
720 fprintf( stderr, "alg = %c\n", alg );
721 ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
722 break;
723 }
724
725 nlen = strlen( tmpres1[0] );
726 if( totallen + nlen > alloclen ) ErrorExit( "LENGTH OVER in Falign\n " );
727 for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
728 for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
729 totallen += nlen;
730 #if 0
731 fprintf( stderr, "%4d\r", totallen );
732 fprintf( stderr, "\n\n" );
733 for( j=0; j<clus1; j++ )
734 {
735 fprintf( stderr, "%s\n", tmpres1[j] );
736 }
737 fprintf( stderr, "-------\n" );
738 for( j=0; j<clus2; j++ )
739 {
740 fprintf( stderr, "%s\n", tmpres2[j] );
741 }
742 #endif
743 }
744 #if KEIKA
745 fprintf( stderr, "DP ... done \n" );
746 #endif
747
748 for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
749 for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
750 #if 0
751 for( j=0; j<clus1; j++ )
752 {
753 fprintf( stderr, "%s\n", result1[j] );
754 }
755 fprintf( stderr, "- - - - - - - - - - -\n" );
756 for( j=0; j<clus2; j++ )
757 {
758 fprintf( stderr, "%s\n", result2[j] );
759 }
760 #endif
761 return( 0.0 );
762 }
763
764
765
Falign(int ** whichmtx,double *** scoringmatrices,double ** n_dynamicmtx,char ** seq1,char ** seq2,double * eff1,double * eff2,double ** eff1s,double ** eff2s,int clus1,int clus2,int alloclen,int * fftlog,int * chudanpt,int chudanref,int * chudanres)766 double Falign( int **whichmtx, double ***scoringmatrices, double **n_dynamicmtx,
767 char **seq1, char **seq2,
768 double *eff1, double *eff2,
769 double **eff1s, double **eff2s,
770 int clus1, int clus2,
771 int alloclen, int *fftlog,
772 int *chudanpt, int chudanref, int *chudanres )
773 {
774 int i, j, k, l, m, maxk;
775 int nlen, nlen2, nlen4;
776 static TLS int prevalloclen = 0;
777 static TLS int crossscoresize = 0;
778 static TLS char **tmpseq1 = NULL;
779 static TLS char **tmpseq2 = NULL;
780 static TLS char **tmpptr1 = NULL;
781 static TLS char **tmpptr2 = NULL;
782 static TLS char **tmpres1 = NULL;
783 static TLS char **tmpres2 = NULL;
784 static TLS char **result1 = NULL;
785 static TLS char **result2 = NULL;
786 #if RND
787 static TLS char **rndseq1 = NULL;
788 static TLS char **rndseq2 = NULL;
789 #endif
790 static TLS Fukusosuu **seqVector1 = NULL;
791 static TLS Fukusosuu **seqVector2 = NULL;
792 static TLS Fukusosuu **naiseki = NULL;
793 static TLS Fukusosuu *naisekiNoWa = NULL;
794 static TLS double *soukan = NULL;
795 static TLS double **crossscore = NULL;
796 int nlentmp;
797 static TLS int *kouho = NULL;
798 static TLS Segment *segment = NULL;
799 static TLS Segment *segment1 = NULL;
800 static TLS Segment *segment2 = NULL;
801 static TLS Segment **sortedseg1 = NULL;
802 static TLS Segment **sortedseg2 = NULL;
803 static TLS int *cut1 = NULL;
804 static TLS int *cut2 = NULL;
805 static TLS char *sgap1, *egap1, *sgap2, *egap2;
806 static TLS int localalloclen = 0;
807 int lag;
808 int tmpint;
809 int count, count0;
810 int len1, len2;
811 int totallen;
812 double totalscore;
813 double dumdb = 0.0;
814 int headgp, tailgp;
815
816
817 if( seq1 == NULL )
818 {
819 if( result1 )
820 {
821 // fprintf( stderr, "Freeing localarrays in Falign\n" );
822 localalloclen = 0;
823 prevalloclen = 0;
824 crossscoresize = 0;
825 mymergesort( 0, 0, NULL );
826 alignableReagion( 0, 0, NULL, NULL, NULL, NULL, NULL );
827 fft( 0, NULL, 1 );
828 A__align( NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 );
829 D__align( NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 );
830 A__align_variousdist( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 );
831 D__align_variousdist( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 );
832 G__align11( NULL, NULL, NULL, 0, 0, 0 );
833 blockAlign2( NULL, NULL, NULL, NULL, NULL, NULL );
834 if( crossscore ) FreeDoubleMtx( crossscore );
835 crossscore = NULL;
836 FreeCharMtx( result1 ); result1 = NULL;
837 FreeCharMtx( result2 );
838 FreeCharMtx( tmpres1 );
839 FreeCharMtx( tmpres2 );
840 FreeCharMtx( tmpseq1 );
841 FreeCharMtx( tmpseq2 );
842 free( sgap1 );
843 free( egap1 );
844 free( sgap2 );
845 free( egap2 );
846 free( kouho );
847 free( cut1 );
848 free( cut2 );
849 free( tmpptr1 );
850 free( tmpptr2 );
851 free( segment );
852 free( segment1 );
853 free( segment2 );
854 free( sortedseg1 );
855 free( sortedseg2 );
856 if( !kobetsubunkatsu )
857 {
858 FreeFukusosuuMtx ( seqVector1 );
859 FreeFukusosuuMtx ( seqVector2 );
860 FreeFukusosuuVec( naisekiNoWa );
861 FreeFukusosuuMtx( naiseki );
862 FreeDoubleVec( soukan );
863 }
864 }
865 else
866 {
867 // fprintf( stderr, "Did not allocate localarrays in Falign\n" );
868 }
869
870 return( 0.0 );
871 }
872
873 len1 = strlen( seq1[0] );
874 len2 = strlen( seq2[0] );
875 nlentmp = MAX( len1, len2 );
876
877 nlen = 1;
878 while( nlentmp >= nlen ) nlen <<= 1;
879 #if 0
880 fprintf( stderr, "### nlen = %d\n", nlen );
881 #endif
882
883 nlen2 = nlen/2; nlen4 = nlen2 / 2;
884
885 #if DEBUG
886 fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
887 fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
888 #endif
889
890 if( prevalloclen != alloclen ) // Falign_noudp mo kaeru
891 {
892 if( prevalloclen )
893 {
894 FreeCharMtx( result1 );
895 FreeCharMtx( result2 );
896 FreeCharMtx( tmpres1 );
897 FreeCharMtx( tmpres2 );
898 }
899 // fprintf( stderr, "\n\n\nreallocating ...\n" );
900 result1 = AllocateCharMtx( njob, alloclen );
901 result2 = AllocateCharMtx( njob, alloclen );
902 tmpres1 = AllocateCharMtx( njob, alloclen );
903 tmpres2 = AllocateCharMtx( njob, alloclen );
904 prevalloclen = alloclen;
905 }
906 if( !localalloclen )
907 {
908 sgap1 = AllocateCharVec( njob );
909 egap1 = AllocateCharVec( njob );
910 sgap2 = AllocateCharVec( njob );
911 egap2 = AllocateCharVec( njob );
912 kouho = AllocateIntVec( NKOUHO );
913 cut1 = AllocateIntVec( MAXSEG );
914 cut2 = AllocateIntVec( MAXSEG );
915 tmpptr1 = AllocateCharMtx( njob, 0 );
916 tmpptr2 = AllocateCharMtx( njob, 0 );
917 // crossscore = AllocateDoubleMtx( MAXSEG, MAXSEG );
918 segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
919 segment1 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
920 segment2 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
921 sortedseg1 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
922 sortedseg2 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
923 if( !( segment && segment1 && segment2 && sortedseg1 && sortedseg2 ) )
924 ErrorExit( "Allocation error\n" );
925
926 if ( scoremtx == -1 ) n20or4or2 = 1;
927 else if( fftscore ) n20or4or2 = 1;
928 else n20or4or2 = 20;
929 }
930 if( localalloclen < nlen )
931 {
932 if( localalloclen )
933 {
934 #if 1
935 if( !kobetsubunkatsu )
936 {
937 FreeFukusosuuMtx ( seqVector1 );
938 FreeFukusosuuMtx ( seqVector2 );
939 FreeFukusosuuVec( naisekiNoWa );
940 FreeFukusosuuMtx( naiseki );
941 FreeDoubleVec( soukan );
942 }
943 FreeCharMtx( tmpseq1 );
944 FreeCharMtx( tmpseq2 );
945 #endif
946 #if RND
947 FreeCharMtx( rndseq1 );
948 FreeCharMtx( rndseq2 );
949 #endif
950 }
951
952 tmpseq1 = AllocateCharMtx( njob, nlen );
953 tmpseq2 = AllocateCharMtx( njob, nlen );
954 if( !kobetsubunkatsu )
955 {
956 naisekiNoWa = AllocateFukusosuuVec( nlen );
957 naiseki = AllocateFukusosuuMtx( n20or4or2, nlen );
958 seqVector1 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
959 seqVector2 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
960 soukan = AllocateDoubleVec( nlen+1 );
961 }
962 #if RND
963 rndseq1 = AllocateCharMtx( njob, nlen );
964 rndseq2 = AllocateCharMtx( njob, nlen );
965 for( i=0; i<njob; i++ )
966 {
967 generateRndSeq( rndseq1[i], nlen );
968 generateRndSeq( rndseq2[i], nlen );
969 }
970 #endif
971 localalloclen = nlen;
972 }
973
974 for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
975 for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );
976
977 #if 0
978 fftfp = fopen( "input_of_Falign", "w" );
979 fprintf( fftfp, "nlen = %d\n", nlen );
980 fprintf( fftfp, "seq1: ( %d sequences ) \n", clus1 );
981 for( i=0; i<clus1; i++ )
982 fprintf( fftfp, "%s\n", seq1[i] );
983 fprintf( fftfp, "seq2: ( %d sequences ) \n", clus2 );
984 for( i=0; i<clus2; i++ )
985 fprintf( fftfp, "%s\n", seq2[i] );
986 fclose( fftfp );
987 system( "less input_of_Falign < /dev/tty > /dev/tty" );
988 #endif
989 if( !kobetsubunkatsu )
990 {
991 if( fftkeika ) fprintf( stderr, " FFT ... " );
992
993 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
994 if( fftscore && scoremtx != -1 )
995 {
996 for( i=0; i<clus1; i++ )
997 {
998 #if 1
999 seq_vec_5( seqVector1[0], polarity, volume, eff1[i], tmpseq1[i] );
1000 #else
1001 seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
1002 seq_vec_2( seqVector1[1], volume, eff1[i], tmpseq1[i] );
1003 #endif
1004 }
1005 }
1006 else
1007 {
1008 #if 0
1009 for( i=0; i<clus1; i++ ) for( j=0; j<n20or4or2; j++ )
1010 seq_vec( seqVector1[j], amino[j], eff1[i], tmpseq1[i] );
1011 #else
1012 for( i=0; i<clus1; i++ )
1013 seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
1014 #endif
1015 }
1016 #if RND
1017 for( i=0; i<clus1; i++ )
1018 {
1019 vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
1020 }
1021 #endif
1022 #if 0
1023 fftfp = fopen( "seqVec", "w" );
1024 fprintf( fftfp, "before transform\n" );
1025 for( k=0; k<n20or4or2; k++ )
1026 {
1027 fprintf( fftfp, "nlen=%d\n", nlen );
1028 fprintf( fftfp, "%c\n", amino[k] );
1029 for( l=0; l<nlen; l++ )
1030 fprintf( fftfp, "%f %f\n", seqVector1[k][l].R, seqVector1[k][l].I );
1031 }
1032 fclose( fftfp );
1033 system( "less seqVec < /dev/tty > /dev/tty" );
1034 #endif
1035
1036 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
1037 if( fftscore && scoremtx != -1 )
1038 {
1039 for( i=0; i<clus2; i++ )
1040 {
1041 #if 1
1042 seq_vec_5( seqVector2[0], polarity, volume, eff2[i], tmpseq2[i] );
1043 #else
1044 seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
1045 seq_vec_2( seqVector2[1], volume, eff2[i], tmpseq2[i] );
1046 #endif
1047 }
1048 }
1049 else
1050 {
1051 #if 0
1052 for( i=0; i<clus2; i++ ) for( j=0; j<n20or4or2; j++ )
1053 seq_vec( seqVector2[j], amino[j], eff2[i], tmpseq2[i] );
1054 #else
1055 for( i=0; i<clus2; i++ )
1056 seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
1057 #endif
1058 }
1059 #if RND
1060 for( i=0; i<clus2; i++ )
1061 {
1062 vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
1063 }
1064 #endif
1065
1066 #if 0
1067 fftfp = fopen( "seqVec2", "w" );
1068 fprintf( fftfp, "before fft\n" );
1069 for( k=0; k<n20or4or2; k++ )
1070 {
1071 fprintf( fftfp, "%c\n", amino[k] );
1072 for( l=0; l<nlen; l++ )
1073 fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
1074 }
1075 fclose( fftfp );
1076 system( "less seqVec2 < /dev/tty > /dev/tty" );
1077 #endif
1078
1079 for( j=0; j<n20or4or2; j++ )
1080 {
1081 fft( nlen, seqVector2[j], 0 );
1082 fft( nlen, seqVector1[j], 0 );
1083 }
1084 #if 0
1085 fftfp = fopen( "seqVec2", "w" );
1086 fprintf( fftfp, "#after fft\n" );
1087 for( k=0; k<n20or4or2; k++ )
1088 {
1089 fprintf( fftfp, "#%c\n", amino[k] );
1090 for( l=0; l<nlen; l++ )
1091 fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
1092 }
1093 fclose( fftfp );
1094 system( "less seqVec2 < /dev/tty > /dev/tty" );
1095 #endif
1096
1097 for( k=0; k<n20or4or2; k++ )
1098 {
1099 for( l=0; l<nlen; l++ )
1100 calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
1101 }
1102 for( l=0; l<nlen; l++ )
1103 {
1104 naisekiNoWa[l].R = 0.0;
1105 naisekiNoWa[l].I = 0.0;
1106 for( k=0; k<n20or4or2; k++ )
1107 {
1108 naisekiNoWa[l].R += naiseki[k][l].R;
1109 naisekiNoWa[l].I += naiseki[k][l].I;
1110 }
1111 }
1112
1113 #if 0
1114 fftfp = fopen( "naisekiNoWa", "w" );
1115 fprintf( fftfp, "#Before fft\n" );
1116 for( l=0; l<nlen; l++ )
1117 fprintf( fftfp, "%d %f %f\n", l, naisekiNoWa[l].R, naisekiNoWa[l].I );
1118 fclose( fftfp );
1119 system( "less naisekiNoWa < /dev/tty > /dev/tty " );
1120 #endif
1121
1122 fft( -nlen, naisekiNoWa, 0 );
1123
1124 for( m=0; m<=nlen2; m++ )
1125 soukan[m] = naisekiNoWa[nlen2-m].R;
1126 for( m=nlen2+1; m<nlen; m++ )
1127 soukan[m] = naisekiNoWa[nlen+nlen2-m].R;
1128
1129 #if 0
1130 fftfp = fopen( "naisekiNoWa", "w" );
1131 fprintf( fftfp, "#After fft\n" );
1132 for( l=0; l<nlen; l++ )
1133 fprintf( fftfp, "%d %f\n", l, naisekiNoWa[l].R );
1134 fclose( fftfp );
1135 fftfp = fopen( "list.plot", "w" );
1136 fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
1137 fclose( fftfp );
1138 system( "/usr/bin/gnuplot list.plot &" );
1139 #endif
1140 #if 0
1141 fprintf( stderr, "soukan\n" );
1142 for( l=0; l<nlen; l++ )
1143 fprintf( stderr, "%d %f\n", l-nlen2, soukan[l] );
1144 #if 0
1145 fftfp = fopen( "list.plot", "w" );
1146 fprintf( fftfp, "plot 'frt'\n pause +1" );
1147 fclose( fftfp );
1148 system( "/usr/bin/gnuplot list.plot" );
1149 #endif
1150 #endif
1151
1152
1153 getKouho( kouho, NKOUHO, soukan, nlen );
1154
1155 #if 0
1156 for( i=0; i<NKOUHO; i++ )
1157 {
1158 fprintf( stderr, "kouho[%d] = %d\n", i, kouho[i] );
1159 }
1160 #endif
1161 }
1162
1163 #if KEIKA
1164 fprintf( stderr, "Searching anchors ... " );
1165 #endif
1166 count = 0;
1167
1168
1169
1170
1171 #define CAND 0
1172 #if CAND
1173 fftfp = fopen( "cand", "w" );
1174 fclose( fftfp );
1175 #endif
1176 if( kobetsubunkatsu )
1177 {
1178 maxk = 1;
1179 kouho[0] = 0;
1180 }
1181 else
1182 {
1183 maxk = NKOUHO;
1184 }
1185
1186 for( k=0; k<maxk; k++ )
1187 {
1188 lag = kouho[k];
1189 if( lag <= -len1 || len2 <= lag ) continue;
1190 zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
1191 #if CAND
1192 fftfp = fopen( "cand", "a" );
1193 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
1194 fprintf( fftfp, "%s\n", tmpptr1[0] );
1195 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
1196 fprintf( fftfp, "%s\n", tmpptr2[0] );
1197 fprintf( fftfp, ">\n", k+1, lag );
1198 fclose( fftfp );
1199 #endif
1200
1201 // fprintf( stderr, "lag = %d\n", lag );
1202 tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );
1203
1204 // if( lag == -50 ) exit( 1 );
1205
1206 if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );
1207
1208
1209 if( tmpint == 0 ) break; // 060430 iinoka ?
1210 while( tmpint-- > 0 )
1211 {
1212 #if 0
1213 if( segment[count].end - segment[count].start < fftWinSize )
1214 {
1215 count++;
1216 continue;
1217 }
1218 #endif
1219 if( lag > 0 )
1220 {
1221 segment1[count].start = segment[count].start ;
1222 segment1[count].end = segment[count].end ;
1223 segment1[count].center = segment[count].center;
1224 segment1[count].score = segment[count].score;
1225
1226 segment2[count].start = segment[count].start + lag;
1227 segment2[count].end = segment[count].end + lag;
1228 segment2[count].center = segment[count].center + lag;
1229 segment2[count].score = segment[count].score ;
1230 }
1231 else
1232 {
1233 segment1[count].start = segment[count].start - lag;
1234 segment1[count].end = segment[count].end - lag;
1235 segment1[count].center = segment[count].center - lag;
1236 segment1[count].score = segment[count].score ;
1237
1238 segment2[count].start = segment[count].start ;
1239 segment2[count].end = segment[count].end ;
1240 segment2[count].center = segment[count].center;
1241 segment2[count].score = segment[count].score ;
1242 }
1243 #if 0
1244 fprintf( stderr, "in 1 %d\n", segment1[count].center );
1245 fprintf( stderr, "in 2 %d\n", segment2[count].center );
1246 #endif
1247 segment1[count].pair = &segment2[count];
1248 segment2[count].pair = &segment1[count];
1249 count++;
1250 }
1251 }
1252 #if 0
1253 if( !kobetsubunkatsu && fftkeika )
1254 fprintf( stderr, "%d anchors found\r", count );
1255 #endif
1256 if( !count && fftNoAnchStop )
1257 ErrorExit( "Cannot detect anchor!" );
1258 #if 0
1259 fprintf( stderr, "RESULT before sort:\n" );
1260 for( l=0; l<count+1; l++ )
1261 {
1262 fprintf( stderr, "cut[%d]=%d, ", l, segment1[l].center );
1263 fprintf( stderr, "%d score = %f\n", segment2[l].center, segment1[l].score );
1264 }
1265 #endif
1266
1267 #if KEIKA
1268 fprintf( stderr, "done. (%d anchors)\n", count );
1269 fprintf( stderr, "Aligning anchors ... " );
1270 #endif
1271 for( i=0; i<count; i++ )
1272 {
1273 sortedseg1[i] = &segment1[i];
1274 sortedseg2[i] = &segment2[i];
1275 }
1276 #if 0
1277 tmpsort( count, sortedseg1 );
1278 tmpsort( count, sortedseg2 );
1279 qsort( sortedseg1, count, sizeof( Segment * ), segcmp );
1280 qsort( sortedseg2, count, sizeof( Segment * ), segcmp );
1281 #else
1282 mymergesort( 0, count-1, sortedseg1 );
1283 mymergesort( 0, count-1, sortedseg2 );
1284 #endif
1285 for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
1286 for( i=0; i<count; i++ ) sortedseg2[i]->number = i;
1287
1288
1289 if( kobetsubunkatsu )
1290 {
1291 for( i=0; i<count; i++ )
1292 {
1293 cut1[i+1] = sortedseg1[i]->center;
1294 cut2[i+1] = sortedseg2[i]->center;
1295 }
1296 cut1[0] = 0;
1297 cut2[0] = 0;
1298 cut1[count+1] = len1;
1299 cut2[count+1] = len2;
1300 count += 2;
1301 }
1302 else
1303 {
1304 if( crossscoresize < count+2 )
1305 {
1306 crossscoresize = count+2;
1307 #if 1
1308 if( fftkeika ) fprintf( stderr, "######allocating crossscore, size = %d\n", crossscoresize );
1309 #endif
1310 if( crossscore ) FreeDoubleMtx( crossscore );
1311 crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
1312 }
1313 for( i=0; i<count+2; i++ ) for( j=0; j<count+2; j++ )
1314 crossscore[i][j] = 0.0;
1315 for( i=0; i<count; i++ )
1316 {
1317 crossscore[segment1[i].number+1][segment1[i].pair->number+1] = segment1[i].score;
1318 cut1[i+1] = sortedseg1[i]->center;
1319 cut2[i+1] = sortedseg2[i]->center;
1320 }
1321
1322 #if 0
1323 fprintf( stderr, "AFTER SORT\n" );
1324 for( i=0; i<count+1; i++ ) fprintf( stderr, "%d, %d\n", cut1[i], cut2[i] );
1325 fprintf( stderr, "crossscore = \n" );
1326 for( i=0; i<count+1; i++ )
1327 {
1328 for( j=0; j<count+1; j++ )
1329 fprintf( stderr, "%.0f ", crossscore[i][j] );
1330 fprintf( stderr, "\n" );
1331 }
1332 #endif
1333
1334 crossscore[0][0] = 10000000.0;
1335 cut1[0] = 0;
1336 cut2[0] = 0;
1337 crossscore[count+1][count+1] = 10000000.0;
1338 cut1[count+1] = len1;
1339 cut2[count+1] = len2;
1340 count += 2;
1341 count0 = count;
1342
1343 blockAlign2( cut1, cut2, sortedseg1, sortedseg2, crossscore, &count );
1344
1345 // if( count-count0 )
1346 // fprintf( stderr, "%d unused anchors\n", count0-count );
1347
1348 if( !kobetsubunkatsu && fftkeika )
1349 fprintf( stderr, "%d anchors found\n", count );
1350 if( fftkeika )
1351 {
1352 if( count0 > count )
1353 {
1354 #if 0
1355 fprintf( stderr, "\7 REPEAT!? \n" );
1356 #else
1357 fprintf( stderr, "REPEAT!? \n" );
1358 #endif
1359 if( fftRepeatStop ) exit( 1 );
1360 }
1361 #if KEIKA
1362 else fprintf( stderr, "done\n" );
1363 #endif
1364 }
1365 }
1366
1367 #if 0
1368 fftfp = fopen( "fft", "a" );
1369 fprintf( fftfp, "RESULT after sort:\n" );
1370 for( l=0; l<count; l++ )
1371 {
1372 fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
1373 fprintf( fftfp, "%d\n", segment2[l].center );
1374 }
1375 fclose( fftfp );
1376 #endif
1377
1378 #if 0
1379 fprintf( stderr, "RESULT after blckalign:\n" );
1380 for( l=0; l<count+1; l++ )
1381 {
1382 fprintf( stderr, "cut : %d %d\n", cut1[l], cut2[l] );
1383 }
1384 #endif
1385
1386 #if 0
1387 fprintf( trap_g, "Devided to %d segments\n", count-1 );
1388 fprintf( trap_g, "%d %d forg\n", MIN( clus1, clus2 ), count-1 );
1389 #endif
1390
1391 totallen = 0;
1392 for( j=0; j<clus1; j++ ) result1[j][0] = 0;
1393 for( j=0; j<clus2; j++ ) result2[j][0] = 0;
1394 totalscore = 0.0;
1395 *fftlog = -1;
1396 for( i=0; i<count-1; i++ )
1397 {
1398 *fftlog += 1;
1399 if( i == 0 ) headgp = outgap; else headgp = 1;
1400 if( i == count-2 ) tailgp = outgap; else tailgp = 1;
1401
1402
1403 if( cut1[i] ) // chuui
1404 {
1405 // getkyokaigap( sgap1, tmpres1, nlen-1, clus1 );
1406 // getkyokaigap( sgap2, tmpres2, nlen-1, clus2 );
1407 getkyokaigap( sgap1, tmpres1, nlen-1, clus1 );
1408 getkyokaigap( sgap2, tmpres2, nlen-1, clus2 );
1409 }
1410 else
1411 {
1412 for( j=0; j<clus1; j++ ) sgap1[j] = 'o';
1413 for( j=0; j<clus2; j++ ) sgap2[j] = 'o';
1414 }
1415 if( cut1[i+1] != len1 ) // chuui
1416 {
1417 getkyokaigap( egap1, seq1, cut1[i+1], clus1 );
1418 getkyokaigap( egap2, seq2, cut2[i+1], clus2 );
1419 }
1420 else
1421 {
1422 for( j=0; j<clus1; j++ ) egap1[j] = 'o';
1423 for( j=0; j<clus2; j++ ) egap2[j] = 'o';
1424 }
1425 #if 0
1426 {
1427 fprintf( stderr, "kyokkaigap1(%d)=", cut1[i]-1 );
1428 for( j=0; j<clus1; j++ )
1429 fprintf( stderr, "%c", sgap1[j] );
1430 fprintf( stderr, "=kyokkaigap1-start\n" );
1431 }
1432 {
1433 fprintf( stderr, "kyokkaigap2(%d)=", cut2[i]-1 );
1434 for( j=0; j<clus2; j++ )
1435 fprintf( stderr, "%c", sgap2[j] );
1436 fprintf( stderr, "=kyokkaigap2-start\n" );
1437 }
1438 {
1439 fprintf( stderr, "kyokkaigap1(%d)=", cut1[i]-1 );
1440 for( j=0; j<clus1; j++ )
1441 fprintf( stderr, "%c", egap1[j] );
1442 fprintf( stderr, "=kyokkaigap1-end\n" );
1443 }
1444 {
1445 fprintf( stderr, "kyokkaigap2(%d)=", cut2[i]-1 );
1446 for( j=0; j<clus2; j++ )
1447 fprintf( stderr, "%c", egap2[j] );
1448 fprintf( stderr, "=kyokkaigap2-end\n" );
1449 }
1450 #endif
1451
1452 #if DEBUG
1453 fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
1454 #else
1455 #if KEIKA
1456 fprintf( stderr, "DP %03d / %03d\r", i+1, count-1 );
1457 #endif
1458 #endif
1459 for( j=0; j<clus1; j++ )
1460 {
1461 strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
1462 tmpres1[j][cut1[i+1]-cut1[i]] = 0;
1463 }
1464 if( kobetsubunkatsu && fftkeika ) commongappick( clus1, tmpres1 ); //dvtditr $B$K8F$P$l$?$H$-(B fftkeika=1
1465 // if( kobetsubunkatsu ) commongappick( clus1, tmpres1 );
1466 for( j=0; j<clus2; j++ )
1467 {
1468 strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
1469 tmpres2[j][cut2[i+1]-cut2[i]] = 0;
1470 }
1471 if( kobetsubunkatsu && fftkeika ) commongappick( clus2, tmpres2 ); //dvtditr $B$K8F$P$l$?$H$-(B fftkeika=1
1472 // if( kobetsubunkatsu ) commongappick( clus2, tmpres2 );
1473
1474 if( constraint )
1475 {
1476 fprintf( stderr, "Not supported\n" );
1477 exit( 1 );
1478 }
1479 #if 0
1480 fprintf( stderr, "i=%d, before alignment", i );
1481 fprintf( stderr, "%4d\n", totallen );
1482 fprintf( stderr, "\n\n" );
1483 for( j=0; j<clus1; j++ )
1484 {
1485 fprintf( stderr, "%s\n", tmpres1[j] );
1486 }
1487 fprintf( stderr, "-------\n" );
1488 for( j=0; j<clus2; j++ )
1489 {
1490 fprintf( stderr, "%s\n", tmpres2[j] );
1491 }
1492 #endif
1493
1494 #if 0
1495 fprintf( stdout, "writing input\n" );
1496 for( j=0; j<clus1; j++ )
1497 {
1498 fprintf( stdout, ">%d of GROUP1\n", j );
1499 fprintf( stdout, "%s\n", tmpres1[j] );
1500 }
1501 for( j=0; j<clus2; j++ )
1502 {
1503 fprintf( stdout, ">%d of GROUP2\n", j );
1504 fprintf( stdout, "%s\n", tmpres2[j] );
1505 }
1506 fflush( stdout );
1507 #endif
1508 switch( alg )
1509 {
1510 case( 'a' ):
1511 totalscore += Aalign( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen );
1512 break;
1513 case( 'M' ):
1514 if( scoringmatrices ) // called by tditeration.c
1515 totalscore += MSalignmm_variousdist( NULL, scoringmatrices, NULL, tmpres1, tmpres2, eff1, eff2, eff1s, eff2s, clus1, clus2, alloclen, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres, headgp, tailgp );
1516 else
1517 totalscore += MSalignmm( n_dynamicmtx, tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres, headgp, tailgp );
1518 // totalscore += MSalignmm( n_dis_consweight_multi, tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres, headgp, tailgp );
1519 break;
1520 case( 'd' ):
1521 if( clus1 == 1 && clus2 == 1 )
1522 {
1523 totalscore += G__align11( n_dynamicmtx, tmpres1, tmpres2, alloclen, headgp, tailgp );
1524 }
1525 else
1526 {
1527 if( scoringmatrices ) // called by tditeration.c
1528 {
1529 totalscore += D__align_variousdist( whichmtx, scoringmatrices, NULL, tmpres1, tmpres2, eff1, eff2, eff1s, eff2s, clus1, clus2, alloclen, NULL, &dumdb, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres, headgp, tailgp );
1530 }
1531 else
1532 totalscore += D__align( n_dynamicmtx, tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumdb, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres, headgp, tailgp );
1533 }
1534 break;
1535 case( 'A' ):
1536 if( clus1 == 1 && clus2 == 1 )
1537 {
1538 totalscore += G__align11( n_dynamicmtx, tmpres1, tmpres2, alloclen, headgp, tailgp );
1539 }
1540 else
1541 {
1542 if( scoringmatrices ) // called by tditeration.c
1543 {
1544 totalscore += A__align_variousdist( whichmtx, scoringmatrices, NULL, tmpres1, tmpres2, eff1, eff2, eff1s, eff2s, clus1, clus2, alloclen, NULL, &dumdb, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres, headgp, tailgp );
1545 }
1546 else
1547 totalscore += A__align( n_dynamicmtx, tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumdb, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres, headgp, tailgp );
1548 }
1549 break;
1550 default:
1551 fprintf( stderr, "alg = %c\n", alg );
1552 ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
1553 break;
1554 }
1555
1556 #ifdef enablemultithread
1557 if( chudanres && *chudanres )
1558 {
1559 // fprintf( stderr, "\n\n## CHUUDAN!!! at Falign_localhom\n" );
1560 return( -1.0 );
1561 }
1562 #endif
1563
1564 nlen = strlen( tmpres1[0] );
1565 if( totallen + nlen > alloclen )
1566 {
1567 fprintf( stderr, "totallen=%d + nlen=%d > alloclen = %d\n", totallen, nlen, alloclen );
1568 ErrorExit( "LENGTH OVER in Falign\n " );
1569 }
1570 for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
1571 for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
1572 totallen += nlen;
1573 #if 0
1574 fprintf( stderr, "$#####$$$$ i=%d", i );
1575 fprintf( stderr, "%4d\n", totallen );
1576 fprintf( stderr, "\n\n" );
1577 for( j=0; j<clus1; j++ )
1578 {
1579 fprintf( stderr, "%s\n", tmpres1[j] );
1580 }
1581 fprintf( stderr, "-------\n" );
1582 for( j=0; j<clus2; j++ )
1583 {
1584 fprintf( stderr, "%s\n", tmpres2[j] );
1585 }
1586 #endif
1587 }
1588 #if KEIKA
1589 fprintf( stderr, "DP ... done \n" );
1590 #endif
1591
1592 for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
1593 for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
1594 #if 0
1595 for( j=0; j<clus1; j++ )
1596 {
1597 fprintf( stderr, "in Falign, %s\n", result1[j] );
1598 }
1599 fprintf( stderr, "- - - - - - - - - - -\n" );
1600 for( j=0; j<clus2; j++ )
1601 {
1602 fprintf( stderr, "in Falign, %s\n", result2[j] );
1603 }
1604 #endif
1605 return( totalscore );
1606 }
1607
1608
1609
1610
1611
1612
1613
1614
1615 /*
1616 sakujo wo kentou (2010/10/05)
1617 */
Falign_udpari_long(int ** whichmtx,double *** scoringmatrices,double ** n_dynamicmtx,char ** seq1,char ** seq2,double * eff1,double * eff2,double ** eff1s,double ** eff2s,int clus1,int clus2,int alloclen,int * fftlog)1618 double Falign_udpari_long(
1619 int **whichmtx, double ***scoringmatrices,
1620 double **n_dynamicmtx,
1621 char **seq1, char **seq2,
1622 double *eff1, double *eff2,
1623 double **eff1s, double **eff2s,
1624 int clus1, int clus2,
1625 int alloclen, int *fftlog )
1626 {
1627 int i, j, k, l, m, maxk;
1628 int nlen, nlen2, nlen4;
1629 static TLS int prevalloclen = 0;
1630 static TLS int crossscoresize = 0;
1631 static TLS char **tmpseq1 = NULL;
1632 static TLS char **tmpseq2 = NULL;
1633 static TLS char **tmpptr1 = NULL;
1634 static TLS char **tmpptr2 = NULL;
1635 static TLS char **tmpres1 = NULL;
1636 static TLS char **tmpres2 = NULL;
1637 static TLS char **result1 = NULL;
1638 static TLS char **result2 = NULL;
1639 #if RND
1640 static TLS char **rndseq1 = NULL;
1641 static TLS char **rndseq2 = NULL;
1642 #endif
1643 static TLS Fukusosuu **seqVector1 = NULL;
1644 static TLS Fukusosuu **seqVector2 = NULL;
1645 static TLS Fukusosuu **naiseki = NULL;
1646 static TLS Fukusosuu *naisekiNoWa = NULL;
1647 static TLS double *soukan = NULL;
1648 static TLS double **crossscore = NULL;
1649 int nlentmp;
1650 static TLS int *kouho = NULL;
1651 static TLS Segment *segment = NULL;
1652 static TLS Segment *segment1 = NULL;
1653 static TLS Segment *segment2 = NULL;
1654 static TLS Segment **sortedseg1 = NULL;
1655 static TLS Segment **sortedseg2 = NULL;
1656 static TLS int *cut1 = NULL;
1657 static TLS int *cut2 = NULL;
1658 static TLS char *sgap1, *egap1, *sgap2, *egap2;
1659 static TLS int localalloclen = 0;
1660 int lag;
1661 int tmpint;
1662 int count, count0;
1663 int len1, len2;
1664 int totallen;
1665 double totalscore;
1666 int nkouho = 0;
1667 int headgp, tailgp;
1668 // double dumfl = 0.0;
1669
1670 if( seq1 == NULL )
1671 {
1672 if( result1 )
1673 {
1674 // fprintf( stderr, "### Freeing localarrays in Falign\n" );
1675 localalloclen = 0;
1676 prevalloclen = 0;
1677 crossscoresize = 0;
1678 mymergesort( 0, 0, NULL );
1679 alignableReagion( 0, 0, NULL, NULL, NULL, NULL, NULL );
1680 fft( 0, NULL, 1 );
1681 A__align( NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 );
1682 A__align_variousdist( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 );
1683 D__align_variousdist( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 );
1684 G__align11( NULL, NULL, NULL, 0, 0, 0 );
1685 blockAlign2( NULL, NULL, NULL, NULL, NULL, NULL );
1686 if( crossscore ) FreeDoubleMtx( crossscore );
1687 crossscore = NULL; // reallocate sareru kanousei ga arunode.
1688 FreeCharMtx( result1 ); result1 = NULL;
1689 FreeCharMtx( result2 );
1690 FreeCharMtx( tmpres1 );
1691 FreeCharMtx( tmpres2 );
1692 FreeCharMtx( tmpseq1 );
1693 FreeCharMtx( tmpseq2 );
1694 free( sgap1 );
1695 free( egap1 );
1696 free( sgap2 );
1697 free( egap2 );
1698 free( kouho );
1699 free( cut1 );
1700 free( cut2 );
1701 free( tmpptr1 );
1702 free( tmpptr2 );
1703 free( segment );
1704 free( segment1 );
1705 free( segment2 );
1706 free( sortedseg1 );
1707 free( sortedseg2 );
1708 if( !kobetsubunkatsu )
1709 {
1710 FreeFukusosuuMtx ( seqVector1 );
1711 FreeFukusosuuMtx ( seqVector2 );
1712 FreeFukusosuuVec( naisekiNoWa );
1713 FreeFukusosuuMtx( naiseki );
1714 FreeDoubleVec( soukan );
1715 }
1716 }
1717 else
1718 {
1719 // fprintf( stderr, "Did not allocate localarrays in Falign\n" );
1720 }
1721
1722 return( 0.0 );
1723 }
1724
1725
1726
1727 len1 = strlen( seq1[0] );
1728 len2 = strlen( seq2[0] );
1729 nlentmp = MAX( len1, len2 );
1730
1731 nlen = 1;
1732 while( nlentmp >= nlen ) nlen <<= 1;
1733 #if 0
1734 fprintf( stderr, "### nlen = %d\n", nlen );
1735 #endif
1736
1737 nlen2 = nlen/2; nlen4 = nlen2 / 2;
1738
1739 #if 0
1740 fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
1741 fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
1742 #endif
1743
1744 if( prevalloclen != alloclen ) // Falign_noudp mo kaeru
1745 {
1746 if( prevalloclen )
1747 {
1748 FreeCharMtx( result1 );
1749 FreeCharMtx( result2 );
1750 FreeCharMtx( tmpres1 );
1751 FreeCharMtx( tmpres2 );
1752 }
1753 // fprintf( stderr, "\n\n\nreallocating ...\n" );
1754 result1 = AllocateCharMtx( njob, alloclen );
1755 result2 = AllocateCharMtx( njob, alloclen );
1756 tmpres1 = AllocateCharMtx( njob, alloclen );
1757 tmpres2 = AllocateCharMtx( njob, alloclen );
1758 prevalloclen = alloclen;
1759 }
1760
1761 if( !localalloclen )
1762 {
1763 sgap1 = AllocateCharVec( njob );
1764 egap1 = AllocateCharVec( njob );
1765 sgap2 = AllocateCharVec( njob );
1766 egap2 = AllocateCharVec( njob );
1767 kouho = AllocateIntVec( NKOUHO_LONG );
1768 cut1 = AllocateIntVec( MAXSEG );
1769 cut2 = AllocateIntVec( MAXSEG );
1770 tmpptr1 = AllocateCharMtx( njob, 0 );
1771 tmpptr2 = AllocateCharMtx( njob, 0 );
1772 segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
1773 segment1 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
1774 segment2 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
1775 sortedseg1 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
1776 sortedseg2 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
1777 if( !( segment && segment1 && segment2 && sortedseg1 && sortedseg2 ) )
1778 ErrorExit( "Allocation error\n" );
1779
1780 if ( scoremtx == -1 ) n20or4or2 = 1;
1781 else if( fftscore ) n20or4or2 = 1;
1782 else n20or4or2 = 20;
1783 }
1784 if( localalloclen < nlen )
1785 {
1786 if( localalloclen )
1787 {
1788 #if 1
1789 if( !kobetsubunkatsu )
1790 {
1791 FreeFukusosuuMtx ( seqVector1 );
1792 FreeFukusosuuMtx ( seqVector2 );
1793 FreeFukusosuuVec( naisekiNoWa );
1794 FreeFukusosuuMtx( naiseki );
1795 FreeDoubleVec( soukan );
1796 }
1797 FreeCharMtx( tmpseq1 );
1798 FreeCharMtx( tmpseq2 );
1799 #endif
1800 #if RND
1801 FreeCharMtx( rndseq1 );
1802 FreeCharMtx( rndseq2 );
1803 #endif
1804 }
1805
1806
1807 tmpseq1 = AllocateCharMtx( njob, nlen );
1808 tmpseq2 = AllocateCharMtx( njob, nlen );
1809 if( !kobetsubunkatsu )
1810 {
1811 naisekiNoWa = AllocateFukusosuuVec( nlen );
1812 naiseki = AllocateFukusosuuMtx( n20or4or2, nlen );
1813 seqVector1 = AllocateFukusosuuMtx( n20or4or2, nlen+1 );
1814 seqVector2 = AllocateFukusosuuMtx( n20or4or2, nlen+1 );
1815 soukan = AllocateDoubleVec( nlen+1 );
1816 }
1817 #if RND
1818 rndseq1 = AllocateCharMtx( njob, nlen );
1819 rndseq2 = AllocateCharMtx( njob, nlen );
1820 for( i=0; i<njob; i++ )
1821 {
1822 generateRndSeq( rndseq1[i], nlen );
1823 generateRndSeq( rndseq2[i], nlen );
1824 }
1825 #endif
1826 localalloclen = nlen;
1827 }
1828
1829 for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
1830 for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );
1831
1832 #if 0
1833 fftfp = fopen( "input_of_Falign", "w" );
1834 fprintf( fftfp, "nlen = %d\n", nlen );
1835 fprintf( fftfp, "seq1: ( %d sequences ) \n", clus1 );
1836 for( i=0; i<clus1; i++ )
1837 fprintf( fftfp, "%s\n", seq1[i] );
1838 fprintf( fftfp, "seq2: ( %d sequences ) \n", clus2 );
1839 for( i=0; i<clus2; i++ )
1840 fprintf( fftfp, "%s\n", seq2[i] );
1841 fclose( fftfp );
1842 system( "less input_of_Falign < /dev/tty > /dev/tty" );
1843 #endif
1844 if( !kobetsubunkatsu )
1845 {
1846 fprintf( stderr, " FFT ... " );
1847
1848 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
1849 if( scoremtx == -1 )
1850 {
1851 for( i=0; i<clus1; i++ )
1852 seq_vec_4( seqVector1[0], eff1[i], tmpseq1[i] );
1853 }
1854 else if( fftscore )
1855 {
1856 for( i=0; i<clus1; i++ )
1857 {
1858 #if 0
1859 seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
1860 seq_vec_2( seqVector1[1], volume, eff1[i], tmpseq1[i] );
1861 #else
1862 seq_vec_5( seqVector1[0], polarity, volume, eff1[i], tmpseq1[i] );
1863 #endif
1864 }
1865 }
1866 else
1867 {
1868 for( i=0; i<clus1; i++ )
1869 seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
1870 }
1871 #if RND
1872 for( i=0; i<clus1; i++ )
1873 {
1874 vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
1875 }
1876 #endif
1877 #if 0
1878 fftfp = fopen( "seqVec", "w" );
1879 fprintf( fftfp, "before transform\n" );
1880 for( k=0; k<n20or4or2; k++ )
1881 {
1882 fprintf( fftfp, "nlen=%d\n", nlen );
1883 fprintf( fftfp, "%c\n", amino[k] );
1884 for( l=0; l<nlen; l++ )
1885 fprintf( fftfp, "%f %f\n", seqVector1[k][l].R, seqVector1[k][l].I );
1886 }
1887 fclose( fftfp );
1888 system( "less seqVec < /dev/tty > /dev/tty" );
1889 #endif
1890
1891 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
1892 if( scoremtx == -1 )
1893 {
1894 for( i=0; i<clus2; i++ )
1895 seq_vec_4( seqVector2[0], eff2[i], tmpseq2[i] );
1896 }
1897 else if( fftscore )
1898 {
1899 for( i=0; i<clus2; i++ )
1900 {
1901 #if 0
1902 seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
1903 seq_vec_2( seqVector2[1], volume, eff2[i], tmpseq2[i] );
1904 #else
1905 seq_vec_5( seqVector2[0], polarity, volume, eff2[i], tmpseq2[i] );
1906 #endif
1907 }
1908 }
1909 else
1910 {
1911 for( i=0; i<clus2; i++ )
1912 seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
1913 }
1914 #if RND
1915 for( i=0; i<clus2; i++ )
1916 {
1917 vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
1918 }
1919 #endif
1920
1921 #if 0
1922 fftfp = fopen( "seqVec2", "w" );
1923 fprintf( fftfp, "before fft\n" );
1924 for( k=0; k<n20or4or2; k++ )
1925 {
1926 fprintf( fftfp, "%c\n", amino[k] );
1927 for( l=0; l<nlen; l++ )
1928 fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
1929 }
1930 fclose( fftfp );
1931 system( "less seqVec2 < /dev/tty > /dev/tty" );
1932 #endif
1933
1934 for( j=0; j<n20or4or2; j++ )
1935 {
1936 fft( nlen, seqVector2[j], 0 );
1937 fft( nlen, seqVector1[j], 0 );
1938 }
1939 #if 0
1940 fftfp = fopen( "seqVec2", "w" );
1941 fprintf( fftfp, "#after fft\n" );
1942 for( k=0; k<n20or4or2; k++ )
1943 {
1944 fprintf( fftfp, "#%c\n", amino[k] );
1945 for( l=0; l<nlen; l++ )
1946 fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
1947 }
1948 fclose( fftfp );
1949 system( "less seqVec2 < /dev/tty > /dev/tty" );
1950 #endif
1951
1952 for( k=0; k<n20or4or2; k++ )
1953 {
1954 for( l=0; l<nlen; l++ )
1955 calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
1956 }
1957 for( l=0; l<nlen; l++ )
1958 {
1959 naisekiNoWa[l].R = 0.0;
1960 naisekiNoWa[l].I = 0.0;
1961 for( k=0; k<n20or4or2; k++ )
1962 {
1963 naisekiNoWa[l].R += naiseki[k][l].R;
1964 naisekiNoWa[l].I += naiseki[k][l].I;
1965 }
1966 }
1967
1968 #if 0
1969 fftfp = fopen( "naisekiNoWa", "w" );
1970 fprintf( fftfp, "#Before fft\n" );
1971 for( l=0; l<nlen; l++ )
1972 fprintf( fftfp, "%d %f %f\n", l, naisekiNoWa[l].R, naisekiNoWa[l].I );
1973 fclose( fftfp );
1974 system( "less naisekiNoWa < /dev/tty > /dev/tty " );
1975 #endif
1976
1977 fft( -nlen, naisekiNoWa, 0 );
1978
1979 for( m=0; m<=nlen2; m++ )
1980 soukan[m] = naisekiNoWa[nlen2-m].R;
1981 for( m=nlen2+1; m<nlen; m++ )
1982 soukan[m] = naisekiNoWa[nlen+nlen2-m].R;
1983
1984 #if 0
1985 fftfp = fopen( "naisekiNoWa", "w" );
1986 fprintf( fftfp, "#After fft\n" );
1987 for( l=0; l<nlen; l++ )
1988 fprintf( fftfp, "%d %f\n", l, naisekiNoWa[l].R );
1989 fclose( fftfp );
1990 fftfp = fopen( "list.plot", "w" );
1991 fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
1992 fclose( fftfp );
1993 system( "/usr/bin/gnuplot list.plot &" );
1994 #endif
1995 #if 0
1996 fprintf( stderr, "soukan\n" );
1997 for( l=0; l<nlen; l++ )
1998 fprintf( stderr, "%d %f\n", l-nlen2, soukan[l] );
1999 #if 0
2000 fftfp = fopen( "list.plot", "w" );
2001 fprintf( fftfp, "plot 'frt'\n pause +1" );
2002 fclose( fftfp );
2003 system( "/usr/bin/gnuplot list.plot" );
2004 #endif
2005 #endif
2006
2007
2008 nkouho = getKouho( kouho, NKOUHO_LONG, soukan, nlen );
2009
2010 #if 0
2011 for( i=0; i<nkouho; i++ )
2012 {
2013 fprintf( stderr, "kouho[%d] = %d\n", i, kouho[i] );
2014 }
2015 #endif
2016 }
2017
2018 #if KEIKA
2019 fprintf( stderr, "Searching anchors ... " );
2020 #endif
2021 count = 0;
2022
2023
2024
2025 #define CAND 0
2026 #if CAND
2027 fftfp = fopen( "cand", "w" );
2028 fclose( fftfp );
2029 #endif
2030 if( kobetsubunkatsu )
2031 {
2032 maxk = 1;
2033 kouho[0] = 0;
2034 }
2035 else
2036 {
2037 maxk = nkouho;
2038 }
2039
2040 for( k=0; k<maxk; k++ )
2041 {
2042 lag = kouho[k];
2043 if( lag <= -len1 || len2 <= lag ) continue;
2044 // fprintf( stderr, "k=%d, lag=%d\n", k, lag );
2045 zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
2046 #if CAND
2047 fftfp = fopen( "cand", "a" );
2048 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
2049 fprintf( fftfp, "%s\n", tmpptr1[0] );
2050 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
2051 fprintf( fftfp, "%s\n", tmpptr2[0] );
2052 fprintf( fftfp, ">\n", k+1, lag );
2053 fclose( fftfp );
2054 #endif
2055
2056 // fprintf( stderr, "lag = %d\n", lag );
2057 tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );
2058 // fprintf( stderr, "lag = %d, %d found\n", lag, tmpint );
2059
2060 // if( lag == -50 ) exit( 1 );
2061
2062 if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );
2063
2064 // fprintf( stderr, "##### k=%d / %d\n", k, maxk );
2065 // if( tmpint == 0 ) break; // 060430 iinoka ? // 090530 yameta
2066 while( tmpint-- > 0 )
2067 {
2068 #if 0
2069 if( segment[count].end - segment[count].start < fftWinSize )
2070 {
2071 count++;
2072 continue;
2073 }
2074 #endif
2075 if( lag > 0 )
2076 {
2077 segment1[count].start = segment[count].start ;
2078 segment1[count].end = segment[count].end ;
2079 segment1[count].center = segment[count].center;
2080 segment1[count].score = segment[count].score;
2081
2082 segment2[count].start = segment[count].start + lag;
2083 segment2[count].end = segment[count].end + lag;
2084 segment2[count].center = segment[count].center + lag;
2085 segment2[count].score = segment[count].score ;
2086 }
2087 else
2088 {
2089 segment1[count].start = segment[count].start - lag;
2090 segment1[count].end = segment[count].end - lag;
2091 segment1[count].center = segment[count].center - lag;
2092 segment1[count].score = segment[count].score ;
2093
2094 segment2[count].start = segment[count].start ;
2095 segment2[count].end = segment[count].end ;
2096 segment2[count].center = segment[count].center;
2097 segment2[count].score = segment[count].score ;
2098 }
2099 #if 0
2100 fprintf( stderr, "##### k=%d / %d\n", k, maxk );
2101 fprintf( stderr, "anchor %d, score = %f\n", count, segment1[count].score );
2102 fprintf( stderr, "in 1 %d\n", segment1[count].center );
2103 fprintf( stderr, "in 2 %d\n", segment2[count].center );
2104 #endif
2105 segment1[count].pair = &segment2[count];
2106 segment2[count].pair = &segment1[count];
2107 count++;
2108 #if 0
2109 fprintf( stderr, "count=%d\n", count );
2110 #endif
2111 }
2112 }
2113 #if 1
2114 if( !kobetsubunkatsu )
2115 fprintf( stderr, "done. (%d anchors) ", count );
2116 #endif
2117 if( !count && fftNoAnchStop )
2118 ErrorExit( "Cannot detect anchor!" );
2119 #if 0
2120 fprintf( stderr, "RESULT before sort:\n" );
2121 for( l=0; l<count+1; l++ )
2122 {
2123 fprintf( stderr, "cut[%d]=%d, ", l, segment1[l].center );
2124 fprintf( stderr, "%d score = %f\n", segment2[l].center, segment1[l].score );
2125 }
2126 #endif
2127
2128 for( i=0; i<count; i++ )
2129 {
2130 sortedseg1[i] = &segment1[i];
2131 sortedseg2[i] = &segment2[i];
2132 }
2133 #if 0
2134 tmpsort( count, sortedseg1 );
2135 tmpsort( count, sortedseg2 );
2136 qsort( sortedseg1, count, sizeof( Segment * ), segcmp );
2137 qsort( sortedseg2, count, sizeof( Segment * ), segcmp );
2138 #else
2139 mymergesort( 0, count-1, sortedseg1 );
2140 mymergesort( 0, count-1, sortedseg2 );
2141 #endif
2142 for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
2143 for( i=0; i<count; i++ ) sortedseg2[i]->number = i;
2144
2145
2146
2147 if( kobetsubunkatsu )
2148 {
2149 for( i=0; i<count; i++ )
2150 {
2151 cut1[i+1] = sortedseg1[i]->center;
2152 cut2[i+1] = sortedseg2[i]->center;
2153 }
2154 cut1[0] = 0;
2155 cut2[0] = 0;
2156 cut1[count+1] = len1;
2157 cut2[count+1] = len2;
2158 count += 2;
2159 }
2160
2161 else
2162 {
2163 if( count < 5000 )
2164 {
2165 if( crossscoresize < count+2 )
2166 {
2167 crossscoresize = count+2;
2168 #if 1
2169 if( fftkeika ) fprintf( stderr, "######allocating crossscore, size = %d\n", crossscoresize );
2170 #endif
2171 if( crossscore ) FreeDoubleMtx( crossscore );
2172 crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
2173 }
2174 for( i=0; i<count+2; i++ ) for( j=0; j<count+2; j++ )
2175 crossscore[i][j] = 0.0;
2176 for( i=0; i<count; i++ )
2177 {
2178 crossscore[segment1[i].number+1][segment1[i].pair->number+1] = segment1[i].score;
2179 cut1[i+1] = sortedseg1[i]->center;
2180 cut2[i+1] = sortedseg2[i]->center;
2181 }
2182
2183 #if 0
2184 fprintf( stderr, "AFTER SORT\n" );
2185 for( i=0; i<count+1; i++ ) fprintf( stderr, "%d, %d\n", cut1[i], cut2[i] );
2186 fprintf( stderr, "crossscore = \n" );
2187 for( i=0; i<count+1; i++ )
2188 {
2189 for( j=0; j<count+1; j++ )
2190 fprintf( stderr, "%.0f ", crossscore[i][j] );
2191 fprintf( stderr, "\n" );
2192 }
2193 #endif
2194
2195 crossscore[0][0] = 10000000.0;
2196 cut1[0] = 0;
2197 cut2[0] = 0;
2198 crossscore[count+1][count+1] = 10000000.0;
2199 cut1[count+1] = len1;
2200 cut2[count+1] = len2;
2201 count += 2;
2202 count0 = count;
2203
2204 // fprintf( stderr, "\n\n\ncalling blockAlign2\n\n\n\n" );
2205 blockAlign2( cut1, cut2, sortedseg1, sortedseg2, crossscore, &count );
2206
2207 // if( count-count0 )
2208 // fprintf( stderr, "%d unused anchors\n", count0-count );
2209
2210 if( !kobetsubunkatsu && fftkeika )
2211 fprintf( stderr, "%d anchors found\n", count );
2212 if( fftkeika )
2213 {
2214 if( count0 > count )
2215 {
2216 #if 0
2217 fprintf( stderr, "\7 REPEAT!? \n" );
2218 #else
2219 fprintf( stderr, "REPEAT!? \n" );
2220 #endif
2221 if( fftRepeatStop ) exit( 1 );
2222 }
2223 #if KEIKA
2224 else fprintf( stderr, "done\n" );
2225 #endif
2226 }
2227 }
2228
2229
2230 else
2231 {
2232 fprintf( stderr, "\nMany anchors were found. The upper-level DP is skipped.\n\n" );
2233
2234 cut1[0] = 0;
2235 cut2[0] = 0;
2236 count0 = 0;
2237 for( i=0; i<count; i++ )
2238 {
2239 // fprintf( stderr, "i=%d, %d-%d ?\n", i, sortedseg1[i]->center, sortedseg1[i]->pair->center );
2240 if( sortedseg1[i]->center > cut1[count0]
2241 && sortedseg1[i]->pair->center > cut2[count0] )
2242 {
2243 count0++;
2244 cut1[count0] = sortedseg1[i]->center;
2245 cut2[count0] = sortedseg1[i]->pair->center;
2246 }
2247 else
2248 {
2249 if( i && sortedseg1[i]->score > sortedseg1[i-1]->score )
2250 {
2251 if( sortedseg1[i]->center > cut1[count0-1]
2252 && sortedseg1[i]->pair->center > cut2[count0-1] )
2253 {
2254 cut1[count0] = sortedseg1[i]->center;
2255 cut2[count0] = sortedseg1[i]->pair->center;
2256 }
2257 else
2258 {
2259 // count0--;
2260 }
2261 }
2262 }
2263 }
2264 // if( count-count0 )
2265 // fprintf( stderr, "%d anchors unused\n", count-count0 );
2266 cut1[count0+1] = len1;
2267 cut2[count0+1] = len2;
2268 count = count0 + 2;
2269 count0 = count;
2270
2271 }
2272 }
2273
2274 // exit( 0 );
2275
2276 #if 0
2277 fftfp = fopen( "fft", "a" );
2278 fprintf( fftfp, "RESULT after sort:\n" );
2279 for( l=0; l<count; l++ )
2280 {
2281 fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
2282 fprintf( fftfp, "%d\n", segment2[l].center );
2283 }
2284 fclose( fftfp );
2285 #endif
2286
2287 #if 0
2288 fprintf( stderr, "RESULT after blckalign:\n" );
2289 for( l=0; l<count+1; l++ )
2290 {
2291 fprintf( stderr, "cut : %d %d\n", cut1[l], cut2[l] );
2292 }
2293 #endif
2294
2295 #if 0
2296 fprintf( trap_g, "Devided to %d segments\n", count-1 );
2297 fprintf( trap_g, "%d %d forg\n", MIN( clus1, clus2 ), count-1 );
2298 #endif
2299
2300 totallen = 0;
2301 for( j=0; j<clus1; j++ ) result1[j][0] = 0;
2302 for( j=0; j<clus2; j++ ) result2[j][0] = 0;
2303 totalscore = 0.0;
2304 *fftlog = -1;
2305 for( i=0; i<count-1; i++ )
2306 {
2307 *fftlog += 1;
2308 if( i == 0 ) headgp = outgap; else headgp = 1;
2309 if( i == count-2 ) tailgp = outgap; else tailgp = 1;
2310
2311 if( cut1[i] )
2312 {
2313 // getkyokaigap( sgap1, seq1, cut1[i]-1, clus1 );
2314 // getkyokaigap( sgap2, seq2, cut2[i]-1, clus2 );
2315 getkyokaigap( sgap1, tmpres1, nlen-1, clus1 );
2316 getkyokaigap( sgap2, tmpres2, nlen-1, clus2 );
2317 }
2318 else
2319 {
2320 for( j=0; j<clus1; j++ ) sgap1[j] = 'o';
2321 for( j=0; j<clus2; j++ ) sgap2[j] = 'o';
2322 }
2323 if( cut1[i+1] != len1 )
2324 {
2325 getkyokaigap( egap1, seq1, cut1[i+1], clus1 );
2326 getkyokaigap( egap2, seq2, cut2[i+1], clus2 );
2327 }
2328 else
2329 {
2330 for( j=0; j<clus1; j++ ) egap1[j] = 'o';
2331 for( j=0; j<clus2; j++ ) egap2[j] = 'o';
2332 }
2333 #if DEBUG
2334 fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
2335 #else
2336 #if 1
2337 fprintf( stderr, "DP %05d / %05d \b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b", i+1, count-1 );
2338 #endif
2339 #endif
2340 for( j=0; j<clus1; j++ )
2341 {
2342 strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
2343 tmpres1[j][cut1[i+1]-cut1[i]] = 0;
2344 }
2345 if( kobetsubunkatsu && fftkeika ) commongappick( clus1, tmpres1 ); //dvtditr $B$K8F$P$l$?$H$-(B fftkeika=1
2346 // if( kobetsubunkatsu ) commongappick( clus1, tmpres1 );
2347 for( j=0; j<clus2; j++ )
2348 {
2349 // fprintf( stderr, "### cut2[i+1]-cut2[i] = %d\n", cut2[i+1]-cut2[i] );
2350 if( cut2[i+1]-cut2[i] <= 0 )
2351 fprintf( stderr, "### cut2[i+1]=%d, cut2[i]=%d\n", cut2[i+1], cut2[i] );
2352 strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
2353 tmpres2[j][cut2[i+1]-cut2[i]] = 0;
2354 }
2355 if( kobetsubunkatsu && fftkeika ) commongappick( clus2, tmpres2 ); //dvtditr $B$K8F$P$l$?$H$-(B fftkeika=1
2356 // if( kobetsubunkatsu ) commongappick( clus2, tmpres2 );
2357
2358 if( constraint )
2359 {
2360 fprintf( stderr, "Not supported\n" );
2361 exit( 1 );
2362 }
2363 #if 0
2364 fprintf( stderr, "i=%d, before alignment", i );
2365 fprintf( stderr, "%4d\n", totallen );
2366 fprintf( stderr, "\n\n" );
2367 for( j=0; j<clus1; j++ )
2368 {
2369 fprintf( stderr, "%s\n", tmpres1[j] );
2370 }
2371 fprintf( stderr, "-------\n" );
2372 for( j=0; j<clus2; j++ )
2373 {
2374 fprintf( stderr, "%s\n", tmpres2[j] );
2375 }
2376 #endif
2377
2378 #if 0
2379 fprintf( stdout, "writing input\n" );
2380 for( j=0; j<clus1; j++ )
2381 {
2382 fprintf( stdout, ">%d of GROUP1\n", j );
2383 fprintf( stdout, "%s\n", tmpres1[j] );
2384 }
2385 for( j=0; j<clus2; j++ )
2386 {
2387 fprintf( stdout, ">%d of GROUP2\n", j );
2388 fprintf( stdout, "%s\n", tmpres2[j] );
2389 }
2390 fflush( stdout );
2391 #endif
2392 switch( alg )
2393 {
2394 case( 'M' ):
2395 if( scoringmatrices ) // called by tditeration.c
2396 totalscore += MSalignmm_variousdist( NULL, scoringmatrices, NULL, tmpres1, tmpres2, eff1, eff2, eff1s, eff2s, clus1, clus2, alloclen, sgap1, sgap2, egap1, egap2, NULL, 0, NULL, headgp, tailgp );
2397 else
2398 totalscore += MSalignmm( n_dynamicmtx, tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, sgap1, sgap2, egap1, egap2, NULL, 0, NULL, headgp, tailgp );
2399 break;
2400 default:
2401 fprintf( stderr, "alg = %c\n", alg );
2402 ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
2403 break;
2404 }
2405
2406 nlen = strlen( tmpres1[0] );
2407 if( totallen + nlen > alloclen )
2408 {
2409 fprintf( stderr, "totallen=%d + nlen=%d > alloclen = %d\n", totallen, nlen, alloclen );
2410 ErrorExit( "LENGTH OVER in Falign\n " );
2411 }
2412 for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
2413 for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
2414 totallen += nlen;
2415 #if 0
2416 fprintf( stderr, "i=%d", i );
2417 fprintf( stderr, "%4d\n", totallen );
2418 fprintf( stderr, "\n\n" );
2419 for( j=0; j<clus1; j++ )
2420 {
2421 fprintf( stderr, "%s\n", tmpres1[j] );
2422 }
2423 fprintf( stderr, "-------\n" );
2424 for( j=0; j<clus2; j++ )
2425 {
2426 fprintf( stderr, "%s\n", tmpres2[j] );
2427 }
2428 #endif
2429 }
2430 #if KEIKA
2431 fprintf( stderr, "DP ... done \n" );
2432 #endif
2433
2434 for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
2435 for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
2436 #if 0
2437 for( j=0; j<clus1; j++ )
2438 {
2439 fprintf( stderr, "%s\n", result1[j] );
2440 }
2441 fprintf( stderr, "- - - - - - - - - - -\n" );
2442 for( j=0; j<clus2; j++ )
2443 {
2444 fprintf( stderr, "%s\n", result2[j] );
2445 }
2446 #endif
2447 return( totalscore );
2448 }
2449