6 static TLS int n20or4or2;
12 #if RND // by D.Mathog
13 static void generateRndSeq( char *seq, int len )
17 *seq++ = (int)( rnd() * n20or4or2 );
24 static void vec_init( Fukusosuu *result, int nlen )
28 result->R = result->I = 0.0;
34 static void vec_init2( Fukusosuu **result, char *seq, double eff, int st, int ed )
37 for( i=st; i<ed; i++ )
38 result[(int)*seq++][i].R += eff;
42 static void seq_vec_2( Fukusosuu *result, double *score, double incr, char *seq )
45 for( ; *seq; result++ )
47 n = amino_n[(int)*seq++];
48 if( n < 20 && n >= 0 ) result->R += incr * score[n];
50 fprintf( stderr, "n=%d, score=%f, inc=%f R=%f\n",n, score[n], incr * score[n], result->R );
55 static void seq_vec_3( Fukusosuu **result, double incr, char *seq )
61 n = amino_n[(int)*seq++];
62 if( n < n20or4or2 && n >= 0 ) result[n][i].R += incr;
66 static void seq_vec_5( Fukusosuu *result, double *score1, double *score2, double incr, char *seq )
69 for( ; *seq; result++ )
71 n = amino_n[(int)*seq++];
72 if( n > 20 ) continue;
73 result->R += incr * score1[n];
74 result->I += incr * score2[n];
76 fprintf( stderr, "n=%d, score=%f, inc=%f R=%f\n",n, score[n], incr * score[n], result->R );
82 static void seq_vec_4( Fukusosuu *result, double incr, char *seq )
85 for( ; *seq; result++ )
100 static void seq_vec( Fukusosuu *result, char query, double incr, char *seq )
107 if( *seq++ == query ) result->R += incr;
110 fprintf( stderr, "i = %d result->R = %f\n", bk-nlen, (result-1)->R );
115 static int checkRepeat( int num, int *cutpos )
122 if( ( tmp = *cutpos++ ) < buf ) return( 1 );
128 static int segcmp( void *ptr1, void *ptr2 )
131 Segment **seg1 = (Segment **)ptr1;
132 Segment **seg2 = (Segment **)ptr2;
134 return( (*seg1)->center - (*seg2)->center );
136 diff = (*seg1)->center - (*seg2)->center;
137 if( diff ) return( diff );
139 diff = (*seg1)->start - (*seg2)->start;
140 if( diff ) return( diff );
142 diff = (*seg1)->end - (*seg2)->end;
143 if( diff ) return( diff );
145 fprintf( stderr, "USE STABLE SORT !!\n" );
153 static void mymergesort( int first, int last, Segment **seg )
156 static TLS int i, j, k, p;
157 static TLS int allo = 0;
158 static TLS Segment **work = NULL;
162 free( work ); work = NULL;
169 if( work ) free( work );
170 work = (Segment **)calloc( allo / 2 + 1, sizeof( Segment *) );
175 middle = ( first + last ) / 2;
176 mymergesort( first, middle, seg );
177 mymergesort( middle+1, last, seg );
179 for( i=first; i<=middle; i++ ) work[p++] = seg[i];
180 i = middle + 1; j = 0; k = first;
181 while( i <= last && j < p )
183 if( work[j]->center <= seg[i]->center )
184 seg[k++] = work[j++];
188 while( j < p ) seg[k++] = work[j++];
193 double Fgetlag( char **seq1, char **seq2,
194 double *eff1, double *eff2,
195 int clus1, int clus2,
199 int nlen, nlen2, nlen4;
200 static TLS int crossscoresize = 0;
201 static TLS char **tmpseq1 = NULL;
202 static TLS char **tmpseq2 = NULL;
203 static TLS char **tmpptr1 = NULL;
204 static TLS char **tmpptr2 = NULL;
205 static TLS char **tmpres1 = NULL;
206 static TLS char **tmpres2 = NULL;
207 static TLS char **result1 = NULL;
208 static TLS char **result2 = NULL;
210 static TLS char **rndseq1 = NULL;
211 static TLS char **rndseq2 = NULL;
213 static TLS Fukusosuu **seqVector1 = NULL;
214 static TLS Fukusosuu **seqVector2 = NULL;
215 static TLS Fukusosuu **naiseki = NULL;
216 static TLS Fukusosuu *naisekiNoWa = NULL;
217 static TLS double *soukan = NULL;
218 static TLS double **crossscore = NULL;
220 static TLS int *kouho = NULL;
221 static TLS Segment *segment = NULL;
222 static TLS Segment *segment1 = NULL;
223 static TLS Segment *segment2 = NULL;
224 static TLS Segment **sortedseg1 = NULL;
225 static TLS Segment **sortedseg2 = NULL;
226 static TLS int *cut1 = NULL;
227 static TLS int *cut2 = NULL;
228 static TLS int localalloclen = 0;
237 len1 = strlen( seq1[0] );
238 len2 = strlen( seq2[0] );
239 nlentmp = MAX( len1, len2 );
242 while( nlentmp >= nlen ) nlen <<= 1;
244 fprintf( stderr, "### nlen = %d\n", nlen );
247 nlen2 = nlen/2; nlen4 = nlen2 / 2;
250 fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
251 fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
256 kouho = AllocateIntVec( NKOUHO );
257 cut1 = AllocateIntVec( MAXSEG );
258 cut2 = AllocateIntVec( MAXSEG );
259 tmpptr1 = AllocateCharMtx( njob, 0 );
260 tmpptr2 = AllocateCharMtx( njob, 0 );
261 result1 = AllocateCharMtx( njob, alloclen );
262 result2 = AllocateCharMtx( njob, alloclen );
263 tmpres1 = AllocateCharMtx( njob, alloclen );
264 tmpres2 = AllocateCharMtx( njob, alloclen );
265 // crossscore = AllocateDoubleMtx( MAXSEG, MAXSEG );
266 segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
267 segment1 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
268 segment2 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
269 sortedseg1 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
270 sortedseg2 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
271 if( !( segment && segment1 && segment2 && sortedseg1 && sortedseg2 ) )
272 ErrorExit( "Allocation error\n" );
274 if ( scoremtx == -1 ) n20or4or2 = 4;
275 else if( fftscore == 1 ) n20or4or2 = 2;
278 if( localalloclen < nlen )
283 FreeFukusosuuMtx ( seqVector1 );
284 FreeFukusosuuMtx ( seqVector2 );
285 FreeFukusosuuVec( naisekiNoWa );
286 FreeFukusosuuMtx( naiseki );
287 FreeDoubleVec( soukan );
288 FreeCharMtx( tmpseq1 );
289 FreeCharMtx( tmpseq2 );
292 FreeCharMtx( rndseq1 );
293 FreeCharMtx( rndseq2 );
298 tmpseq1 = AllocateCharMtx( njob, nlen );
299 tmpseq2 = AllocateCharMtx( njob, nlen );
300 naisekiNoWa = AllocateFukusosuuVec( nlen );
301 naiseki = AllocateFukusosuuMtx( n20or4or2, nlen );
302 seqVector1 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
303 seqVector2 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
304 soukan = AllocateDoubleVec( nlen+1 );
307 rndseq1 = AllocateCharMtx( njob, nlen );
308 rndseq2 = AllocateCharMtx( njob, nlen );
309 for( i=0; i<njob; i++ )
311 generateRndSeq( rndseq1[i], nlen );
312 generateRndSeq( rndseq2[i], nlen );
315 localalloclen = nlen;
318 for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
319 for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );
322 fftfp = fopen( "input_of_Falign", "w" );
323 fprintf( fftfp, "nlen = %d\n", nlen );
324 fprintf( fftfp, "seq1: ( %d sequences ) \n", clus1 );
325 for( i=0; i<clus1; i++ )
326 fprintf( fftfp, "%s\n", seq1[i] );
327 fprintf( fftfp, "seq2: ( %d sequences ) \n", clus2 );
328 for( i=0; i<clus2; i++ )
329 fprintf( fftfp, "%s\n", seq2[i] );
331 system( "less input_of_Falign < /dev/tty > /dev/tty" );
334 if( fftkeika ) fprintf( stderr, " FFT ... " );
336 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
337 if( fftscore && scoremtx != -1 )
339 for( i=0; i<clus1; i++ )
341 seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
342 seq_vec_2( seqVector1[1], volume, eff1[i], tmpseq1[i] );
348 for( i=0; i<clus1; i++ ) for( j=0; j<n20or4or2; j++ )
349 seq_vec( seqVector1[j], amino[j], eff1[i], tmpseq1[i] );
351 for( i=0; i<clus1; i++ )
352 seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
356 for( i=0; i<clus1; i++ )
358 vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
362 fftfp = fopen( "seqVec", "w" );
363 fprintf( fftfp, "before transform\n" );
364 for( k=0; k<n20or4or2; k++ )
366 fprintf( fftfp, "nlen=%d\n", nlen );
367 fprintf( fftfp, "%c\n", amino[k] );
368 for( l=0; l<nlen; l++ )
369 fprintf( fftfp, "%f %f\n", seqVector1[k][l].R, seqVector1[k][l].I );
372 system( "less seqVec < /dev/tty > /dev/tty" );
375 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
376 if( fftscore && scoremtx != -1 )
378 for( i=0; i<clus2; i++ )
380 seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
381 seq_vec_2( seqVector2[1], volume, eff2[i], tmpseq2[i] );
387 for( i=0; i<clus2; i++ ) for( j=0; j<n20or4or2; j++ )
388 seq_vec( seqVector2[j], amino[j], eff2[i], tmpseq2[i] );
390 for( i=0; i<clus2; i++ )
391 seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
395 for( i=0; i<clus2; i++ )
397 vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
402 fftfp = fopen( "seqVec2", "w" );
403 fprintf( fftfp, "before fft\n" );
404 for( k=0; k<n20or4or2; k++ )
406 fprintf( fftfp, "%c\n", amino[k] );
407 for( l=0; l<nlen; l++ )
408 fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
411 system( "less seqVec2 < /dev/tty > /dev/tty" );
414 for( j=0; j<n20or4or2; j++ )
416 fft( nlen, seqVector2[j], 0 );
417 fft( nlen, seqVector1[j], 0 );
420 fftfp = fopen( "seqVec2", "w" );
421 fprintf( fftfp, "#after fft\n" );
422 for( k=0; k<n20or4or2; k++ )
424 fprintf( fftfp, "#%c\n", amino[k] );
425 for( l=0; l<nlen; l++ )
426 fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
429 system( "less seqVec2 < /dev/tty > /dev/tty" );
432 for( k=0; k<n20or4or2; k++ )
434 for( l=0; l<nlen; l++ )
435 calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
437 for( l=0; l<nlen; l++ )
439 naisekiNoWa[l].R = 0.0;
440 naisekiNoWa[l].I = 0.0;
441 for( k=0; k<n20or4or2; k++ )
443 naisekiNoWa[l].R += naiseki[k][l].R;
444 naisekiNoWa[l].I += naiseki[k][l].I;
449 fftfp = fopen( "naisekiNoWa", "w" );
450 fprintf( fftfp, "#Before fft\n" );
451 for( l=0; l<nlen; l++ )
452 fprintf( fftfp, "%d %f %f\n", l, naisekiNoWa[l].R, naisekiNoWa[l].I );
454 system( "less naisekiNoWa < /dev/tty > /dev/tty " );
457 fft( -nlen, naisekiNoWa, 0 );
459 for( m=0; m<=nlen2; m++ )
460 soukan[m] = naisekiNoWa[nlen2-m].R;
461 for( m=nlen2+1; m<nlen; m++ )
462 soukan[m] = naisekiNoWa[nlen+nlen2-m].R;
465 fftfp = fopen( "naisekiNoWa", "w" );
466 fprintf( fftfp, "#After fft\n" );
467 for( l=0; l<nlen; l++ )
468 fprintf( fftfp, "%d %f\n", l, naisekiNoWa[l].R );
470 fftfp = fopen( "list.plot", "w" );
471 fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
473 system( "/usr/bin/gnuplot list.plot &" );
476 fprintf( stderr, "frt write start\n" );
477 fftfp = fopen( "frt", "w" );
478 for( l=0; l<nlen; l++ )
479 fprintf( fftfp, "%d %f\n", l-nlen2, soukan[l] );
481 system( "less frt < /dev/tty > /dev/tty" );
483 fftfp = fopen( "list.plot", "w" );
484 fprintf( fftfp, "plot 'frt'\n pause +1" );
486 system( "/usr/bin/gnuplot list.plot" );
491 getKouho( kouho, NKOUHO, soukan, nlen );
494 for( i=0; i<NKOUHO; i++ )
496 fprintf( stdout, "kouho[%d] = %d\n", i, kouho[i] );
501 fprintf( stderr, "Searching anchors ... " );
509 fftfp = fopen( "cand", "w" );
513 for( k=0; k<NKOUHO; k++ )
517 zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
519 fftfp = fopen( "cand", "a" );
520 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
521 fprintf( fftfp, "%s\n", tmpptr1[0] );
522 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
523 fprintf( fftfp, "%s\n", tmpptr2[0] );
524 fprintf( fftfp, ">\n", k+1, lag );
527 tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );
529 if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );
532 if( tmpint == 0 ) break; // 060430 iinoka ?
533 while( tmpint-- > 0 )
537 segment1[count].start = segment[count].start ;
538 segment1[count].end = segment[count].end ;
539 segment1[count].center = segment[count].center;
540 segment1[count].score = segment[count].score;
542 segment2[count].start = segment[count].start + lag;
543 segment2[count].end = segment[count].end + lag;
544 segment2[count].center = segment[count].center + lag;
545 segment2[count].score = segment[count].score ;
549 segment1[count].start = segment[count].start - lag;
550 segment1[count].end = segment[count].end - lag;
551 segment1[count].center = segment[count].center - lag;
552 segment1[count].score = segment[count].score ;
554 segment2[count].start = segment[count].start ;
555 segment2[count].end = segment[count].end ;
556 segment2[count].center = segment[count].center;
557 segment2[count].score = segment[count].score ;
560 fprintf( stderr, "Goukaku=%dko\n", tmpint );
561 fprintf( stderr, "in 1 %d\n", segment1[count].center );
562 fprintf( stderr, "in 2 %d\n", segment2[count].center );
564 segment1[count].pair = &segment2[count];
565 segment2[count].pair = &segment1[count];
568 fprintf( stderr, "count=%d\n", count );
574 fprintf( stderr, "done. (%d anchors)\r", count );
576 if( !count && fftNoAnchStop )
577 ErrorExit( "Cannot detect anchor!" );
579 fprintf( stdout, "RESULT before sort:\n" );
580 for( l=0; l<count+1; l++ )
582 fprintf( stdout, "cut[%d]=%d, ", l, segment1[l].center );
583 fprintf( stdout, "%d score = %f\n", segment2[l].center, segment1[l].score );
589 fprintf( stderr, "Aligning anchors ... " );
591 for( i=0; i<count; i++ )
593 sortedseg1[i] = &segment1[i];
594 sortedseg2[i] = &segment2[i];
598 mymergesort( 0, count-1, sortedseg1 );
599 mymergesort( 0, count-1, sortedseg2 );
600 for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
601 for( i=0; i<count; i++ ) sortedseg2[i]->number = i;
603 if( crossscoresize < count+2 )
605 crossscoresize = count+2;
606 fprintf( stderr, "####################################################################################################################################allocating crossscore, size = %d\n", crossscoresize );
607 if( crossscore ) FreeDoubleMtx( crossscore );
608 crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
611 for( i=0; i<count+2; i++ ) for( j=0; j<count+2; j++ )
612 crossscore[i][j] = 0.0;
613 for( i=0; i<count; i++ )
615 crossscore[segment1[i].number+1][segment1[i].pair->number+1] = segment1[i].score;
616 cut1[i+1] = sortedseg1[i]->center;
617 cut2[i+1] = sortedseg2[i]->center;
621 fprintf( stderr, "AFTER SORT\n" );
622 for( i=0; i<count; i++ ) fprintf( stderr, "%d, %d\n", segment1[i].start, segment2[i].start );
625 crossscore[0][0] = 10000000.0;
628 crossscore[count+1][count+1] = 10000000.0;
629 cut1[count+1] = len1;
630 cut2[count+1] = len2;
634 blockAlign2( cut1, cut2, sortedseg1, sortedseg2, crossscore, &count );
640 fprintf( stderr, "REPEAT!? \n" );
641 if( fftRepeatStop ) exit( 1 );
645 fprintf( stderr, "done\n" );
646 fprintf( stderr, "done. (%d anchors)\n", count );
651 fftfp = fopen( "fft", "a" );
652 fprintf( fftfp, "RESULT after sort:\n" );
653 for( l=0; l<count; l++ )
655 fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
656 fprintf( fftfp, "%d\n", segment2[l].center );
662 fftfp = fopen( "fft", "a" );
663 fprintf( fftfp, "RESULT after sort:\n" );
664 for( l=0; l<count; l++ )
666 fprintf( fftfp, "cut : %d %d\n", cut1[l], cut2[l] );
672 fprintf( trap_g, "Devided to %d segments\n", count-1 );
673 fprintf( trap_g, "%d %d forg\n", MIN( clus1, clus2 ), count-1 );
677 for( j=0; j<clus1; j++ ) result1[j][0] = 0;
678 for( j=0; j<clus2; j++ ) result2[j][0] = 0;
679 for( i=0; i<count-1; i++ )
681 if( i == 0 ) headgp = outgap; else headgp = 1;
682 if( i == count-2 ) tailgp = outgap; else tailgp = 1;
685 fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
688 fprintf( stderr, "DP %03d / %03d\r", i+1, count-1 );
691 for( j=0; j<clus1; j++ )
693 strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
694 tmpres1[j][cut1[i+1]-cut1[i]] = 0;
696 for( j=0; j<clus2; j++ )
698 strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
699 tmpres2[j][cut2[i+1]-cut2[i]] = 0;
704 Aalign( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen );
707 MSalignmm( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, headgp, tailgp );
710 if( clus1 == 1 && clus2 == 1 )
711 G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
713 A__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, headgp, tailgp );
716 if( clus1 == 1 && clus2 == 1 )
717 G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
719 H__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
722 if( clus1 == 1 && clus2 == 1 )
723 G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
725 Q__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
728 fprintf( stderr, "alg = %c\n", alg );
729 ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
733 nlen = strlen( tmpres1[0] );
734 if( totallen + nlen > alloclen ) ErrorExit( "LENGTH OVER in Falign\n " );
735 for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
736 for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
739 fprintf( stderr, "%4d\r", totallen );
740 fprintf( stderr, "\n\n" );
741 for( j=0; j<clus1; j++ )
743 fprintf( stderr, "%s\n", tmpres1[j] );
745 fprintf( stderr, "-------\n" );
746 for( j=0; j<clus2; j++ )
748 fprintf( stderr, "%s\n", tmpres2[j] );
753 fprintf( stderr, "DP ... done \n" );
756 for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
757 for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
759 for( j=0; j<clus1; j++ )
761 fprintf( stderr, "%s\n", result1[j] );
763 fprintf( stderr, "- - - - - - - - - - -\n" );
764 for( j=0; j<clus2; j++ )
766 fprintf( stderr, "%s\n", result2[j] );
774 float Falign( char **seq1, char **seq2,
775 double *eff1, double *eff2,
776 int clus1, int clus2,
777 int alloclen, int *fftlog,
778 int *chudanpt, int chudanref, int *chudanres )
780 int i, j, k, l, m, maxk;
781 int nlen, nlen2, nlen4;
782 static TLS int prevalloclen = 0;
783 static TLS int crossscoresize = 0;
784 static TLS char **tmpseq1 = NULL;
785 static TLS char **tmpseq2 = NULL;
786 static TLS char **tmpptr1 = NULL;
787 static TLS char **tmpptr2 = NULL;
788 static TLS char **tmpres1 = NULL;
789 static TLS char **tmpres2 = NULL;
790 static TLS char **result1 = NULL;
791 static TLS char **result2 = NULL;
793 static TLS char **rndseq1 = NULL;
794 static TLS char **rndseq2 = NULL;
796 static TLS Fukusosuu **seqVector1 = NULL;
797 static TLS Fukusosuu **seqVector2 = NULL;
798 static TLS Fukusosuu **naiseki = NULL;
799 static TLS Fukusosuu *naisekiNoWa = NULL;
800 static TLS double *soukan = NULL;
801 static TLS double **crossscore = NULL;
803 static TLS int *kouho = NULL;
804 static TLS Segment *segment = NULL;
805 static TLS Segment *segment1 = NULL;
806 static TLS Segment *segment2 = NULL;
807 static TLS Segment **sortedseg1 = NULL;
808 static TLS Segment **sortedseg2 = NULL;
809 static TLS int *cut1 = NULL;
810 static TLS int *cut2 = NULL;
811 static TLS char *sgap1, *egap1, *sgap2, *egap2;
812 static TLS int localalloclen = 0;
827 // fprintf( stderr, "Freeing localarrays in Falign\n" );
829 mymergesort( 0, 0, NULL );
830 alignableReagion( 0, 0, NULL, NULL, NULL, NULL, NULL );
832 A__align( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 );
833 G__align11( NULL, NULL, 0, 0, 0 );
834 blockAlign2( NULL, NULL, NULL, NULL, NULL, NULL );
835 if( crossscore ) FreeDoubleMtx( crossscore );
836 FreeCharMtx( result1 );
837 FreeCharMtx( result2 );
838 FreeCharMtx( tmpres1 );
839 FreeCharMtx( tmpres2 );
840 FreeCharMtx( tmpseq1 );
841 FreeCharMtx( tmpseq2 );
856 if( !kobetsubunkatsu )
858 FreeFukusosuuMtx ( seqVector1 );
859 FreeFukusosuuMtx ( seqVector2 );
860 FreeFukusosuuVec( naisekiNoWa );
861 FreeFukusosuuMtx( naiseki );
862 FreeDoubleVec( soukan );
867 // fprintf( stderr, "Did not allocate localarrays in Falign\n" );
873 len1 = strlen( seq1[0] );
874 len2 = strlen( seq2[0] );
875 nlentmp = MAX( len1, len2 );
878 while( nlentmp >= nlen ) nlen <<= 1;
880 fprintf( stderr, "### nlen = %d\n", nlen );
883 nlen2 = nlen/2; nlen4 = nlen2 / 2;
886 fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
887 fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
890 if( prevalloclen != alloclen ) // Falign_noudp mo kaeru
894 FreeCharMtx( result1 );
895 FreeCharMtx( result2 );
896 FreeCharMtx( tmpres1 );
897 FreeCharMtx( tmpres2 );
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;
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" );
926 if ( scoremtx == -1 ) n20or4or2 = 1;
927 else if( fftscore ) n20or4or2 = 1;
930 if( localalloclen < nlen )
935 if( !kobetsubunkatsu )
937 FreeFukusosuuMtx ( seqVector1 );
938 FreeFukusosuuMtx ( seqVector2 );
939 FreeFukusosuuVec( naisekiNoWa );
940 FreeFukusosuuMtx( naiseki );
941 FreeDoubleVec( soukan );
943 FreeCharMtx( tmpseq1 );
944 FreeCharMtx( tmpseq2 );
947 FreeCharMtx( rndseq1 );
948 FreeCharMtx( rndseq2 );
952 tmpseq1 = AllocateCharMtx( njob, nlen );
953 tmpseq2 = AllocateCharMtx( njob, nlen );
954 if( !kobetsubunkatsu )
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 );
963 rndseq1 = AllocateCharMtx( njob, nlen );
964 rndseq2 = AllocateCharMtx( njob, nlen );
965 for( i=0; i<njob; i++ )
967 generateRndSeq( rndseq1[i], nlen );
968 generateRndSeq( rndseq2[i], nlen );
971 localalloclen = nlen;
974 for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
975 for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );
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] );
987 system( "less input_of_Falign < /dev/tty > /dev/tty" );
989 if( !kobetsubunkatsu )
991 if( fftkeika ) fprintf( stderr, " FFT ... " );
993 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
994 if( fftscore && scoremtx != -1 )
996 for( i=0; i<clus1; i++ )
999 seq_vec_5( seqVector1[0], polarity, volume, eff1[i], tmpseq1[i] );
1001 seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
1002 seq_vec_2( seqVector1[1], volume, eff1[i], tmpseq1[i] );
1009 for( i=0; i<clus1; i++ ) for( j=0; j<n20or4or2; j++ )
1010 seq_vec( seqVector1[j], amino[j], eff1[i], tmpseq1[i] );
1012 for( i=0; i<clus1; i++ )
1013 seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
1017 for( i=0; i<clus1; i++ )
1019 vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
1023 fftfp = fopen( "seqVec", "w" );
1024 fprintf( fftfp, "before transform\n" );
1025 for( k=0; k<n20or4or2; k++ )
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 );
1033 system( "less seqVec < /dev/tty > /dev/tty" );
1036 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
1037 if( fftscore && scoremtx != -1 )
1039 for( i=0; i<clus2; i++ )
1042 seq_vec_5( seqVector2[0], polarity, volume, eff2[i], tmpseq2[i] );
1044 seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
1045 seq_vec_2( seqVector2[1], volume, eff2[i], tmpseq2[i] );
1052 for( i=0; i<clus2; i++ ) for( j=0; j<n20or4or2; j++ )
1053 seq_vec( seqVector2[j], amino[j], eff2[i], tmpseq2[i] );
1055 for( i=0; i<clus2; i++ )
1056 seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
1060 for( i=0; i<clus2; i++ )
1062 vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
1067 fftfp = fopen( "seqVec2", "w" );
1068 fprintf( fftfp, "before fft\n" );
1069 for( k=0; k<n20or4or2; k++ )
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 );
1076 system( "less seqVec2 < /dev/tty > /dev/tty" );
1079 for( j=0; j<n20or4or2; j++ )
1081 fft( nlen, seqVector2[j], 0 );
1082 fft( nlen, seqVector1[j], 0 );
1085 fftfp = fopen( "seqVec2", "w" );
1086 fprintf( fftfp, "#after fft\n" );
1087 for( k=0; k<n20or4or2; k++ )
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 );
1094 system( "less seqVec2 < /dev/tty > /dev/tty" );
1097 for( k=0; k<n20or4or2; k++ )
1099 for( l=0; l<nlen; l++ )
1100 calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
1102 for( l=0; l<nlen; l++ )
1104 naisekiNoWa[l].R = 0.0;
1105 naisekiNoWa[l].I = 0.0;
1106 for( k=0; k<n20or4or2; k++ )
1108 naisekiNoWa[l].R += naiseki[k][l].R;
1109 naisekiNoWa[l].I += naiseki[k][l].I;
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 );
1119 system( "less naisekiNoWa < /dev/tty > /dev/tty " );
1122 fft( -nlen, naisekiNoWa, 0 );
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;
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 );
1135 fftfp = fopen( "list.plot", "w" );
1136 fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
1138 system( "/usr/bin/gnuplot list.plot &" );
1141 fprintf( stderr, "soukan\n" );
1142 for( l=0; l<nlen; l++ )
1143 fprintf( stderr, "%d %f\n", l-nlen2, soukan[l] );
1145 fftfp = fopen( "list.plot", "w" );
1146 fprintf( fftfp, "plot 'frt'\n pause +1" );
1148 system( "/usr/bin/gnuplot list.plot" );
1153 getKouho( kouho, NKOUHO, soukan, nlen );
1156 for( i=0; i<NKOUHO; i++ )
1158 fprintf( stderr, "kouho[%d] = %d\n", i, kouho[i] );
1164 fprintf( stderr, "Searching anchors ... " );
1172 fftfp = fopen( "cand", "w" );
1175 if( kobetsubunkatsu )
1185 for( k=0; k<maxk; k++ )
1188 if( lag <= -len1 || len2 <= lag ) continue;
1189 zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
1191 fftfp = fopen( "cand", "a" );
1192 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
1193 fprintf( fftfp, "%s\n", tmpptr1[0] );
1194 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
1195 fprintf( fftfp, "%s\n", tmpptr2[0] );
1196 fprintf( fftfp, ">\n", k+1, lag );
1200 // fprintf( stderr, "lag = %d\n", lag );
1201 tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );
1203 // if( lag == -50 ) exit( 1 );
1205 if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );
1208 if( tmpint == 0 ) break; // 060430 iinoka ?
1209 while( tmpint-- > 0 )
1212 if( segment[count].end - segment[count].start < fftWinSize )
1220 segment1[count].start = segment[count].start ;
1221 segment1[count].end = segment[count].end ;
1222 segment1[count].center = segment[count].center;
1223 segment1[count].score = segment[count].score;
1225 segment2[count].start = segment[count].start + lag;
1226 segment2[count].end = segment[count].end + lag;
1227 segment2[count].center = segment[count].center + lag;
1228 segment2[count].score = segment[count].score ;
1232 segment1[count].start = segment[count].start - lag;
1233 segment1[count].end = segment[count].end - lag;
1234 segment1[count].center = segment[count].center - lag;
1235 segment1[count].score = segment[count].score ;
1237 segment2[count].start = segment[count].start ;
1238 segment2[count].end = segment[count].end ;
1239 segment2[count].center = segment[count].center;
1240 segment2[count].score = segment[count].score ;
1243 fprintf( stderr, "in 1 %d\n", segment1[count].center );
1244 fprintf( stderr, "in 2 %d\n", segment2[count].center );
1246 segment1[count].pair = &segment2[count];
1247 segment2[count].pair = &segment1[count];
1252 if( !kobetsubunkatsu && fftkeika )
1253 fprintf( stderr, "%d anchors found\r", count );
1255 if( !count && fftNoAnchStop )
1256 ErrorExit( "Cannot detect anchor!" );
1258 fprintf( stderr, "RESULT before sort:\n" );
1259 for( l=0; l<count+1; l++ )
1261 fprintf( stderr, "cut[%d]=%d, ", l, segment1[l].center );
1262 fprintf( stderr, "%d score = %f\n", segment2[l].center, segment1[l].score );
1267 fprintf( stderr, "done. (%d anchors)\n", count );
1268 fprintf( stderr, "Aligning anchors ... " );
1270 for( i=0; i<count; i++ )
1272 sortedseg1[i] = &segment1[i];
1273 sortedseg2[i] = &segment2[i];
1276 tmpsort( count, sortedseg1 );
1277 tmpsort( count, sortedseg2 );
1278 qsort( sortedseg1, count, sizeof( Segment * ), segcmp );
1279 qsort( sortedseg2, count, sizeof( Segment * ), segcmp );
1281 mymergesort( 0, count-1, sortedseg1 );
1282 mymergesort( 0, count-1, sortedseg2 );
1284 for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
1285 for( i=0; i<count; i++ ) sortedseg2[i]->number = i;
1288 if( kobetsubunkatsu )
1290 for( i=0; i<count; i++ )
1292 cut1[i+1] = sortedseg1[i]->center;
1293 cut2[i+1] = sortedseg2[i]->center;
1297 cut1[count+1] = len1;
1298 cut2[count+1] = len2;
1303 if( crossscoresize < count+2 )
1305 crossscoresize = count+2;
1307 if( fftkeika ) fprintf( stderr, "######allocating crossscore, size = %d\n", crossscoresize );
1309 if( crossscore ) FreeDoubleMtx( crossscore );
1310 crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
1312 for( i=0; i<count+2; i++ ) for( j=0; j<count+2; j++ )
1313 crossscore[i][j] = 0.0;
1314 for( i=0; i<count; i++ )
1316 crossscore[segment1[i].number+1][segment1[i].pair->number+1] = segment1[i].score;
1317 cut1[i+1] = sortedseg1[i]->center;
1318 cut2[i+1] = sortedseg2[i]->center;
1322 fprintf( stderr, "AFTER SORT\n" );
1323 for( i=0; i<count+1; i++ ) fprintf( stderr, "%d, %d\n", cut1[i], cut2[i] );
1324 fprintf( stderr, "crossscore = \n" );
1325 for( i=0; i<count+1; i++ )
1327 for( j=0; j<count+1; j++ )
1328 fprintf( stderr, "%.0f ", crossscore[i][j] );
1329 fprintf( stderr, "\n" );
1333 crossscore[0][0] = 10000000.0;
1336 crossscore[count+1][count+1] = 10000000.0;
1337 cut1[count+1] = len1;
1338 cut2[count+1] = len2;
1342 blockAlign2( cut1, cut2, sortedseg1, sortedseg2, crossscore, &count );
1344 // if( count-count0 )
1345 // fprintf( stderr, "%d unused anchors\n", count0-count );
1347 if( !kobetsubunkatsu && fftkeika )
1348 fprintf( stderr, "%d anchors found\n", count );
1351 if( count0 > count )
1354 fprintf( stderr, "\7 REPEAT!? \n" );
1356 fprintf( stderr, "REPEAT!? \n" );
1358 if( fftRepeatStop ) exit( 1 );
1361 else fprintf( stderr, "done\n" );
1367 fftfp = fopen( "fft", "a" );
1368 fprintf( fftfp, "RESULT after sort:\n" );
1369 for( l=0; l<count; l++ )
1371 fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
1372 fprintf( fftfp, "%d\n", segment2[l].center );
1378 fprintf( stderr, "RESULT after blckalign:\n" );
1379 for( l=0; l<count+1; l++ )
1381 fprintf( stderr, "cut : %d %d\n", cut1[l], cut2[l] );
1386 fprintf( trap_g, "Devided to %d segments\n", count-1 );
1387 fprintf( trap_g, "%d %d forg\n", MIN( clus1, clus2 ), count-1 );
1391 for( j=0; j<clus1; j++ ) result1[j][0] = 0;
1392 for( j=0; j<clus2; j++ ) result2[j][0] = 0;
1395 for( i=0; i<count-1; i++ )
1398 if( i == 0 ) headgp = outgap; else headgp = 1;
1399 if( i == count-2 ) tailgp = outgap; else tailgp = 1;
1402 if( cut1[i] ) // chuui
1404 // getkyokaigap( sgap1, tmpres1, nlen-1, clus1 );
1405 // getkyokaigap( sgap2, tmpres2, nlen-1, clus2 );
1406 getkyokaigap( sgap1, tmpres1, nlen-1, clus1 );
1407 getkyokaigap( sgap2, tmpres2, nlen-1, clus2 );
1411 for( j=0; j<clus1; j++ ) sgap1[j] = 'o';
1412 for( j=0; j<clus2; j++ ) sgap2[j] = 'o';
1414 if( cut1[i+1] != len1 ) // chuui
1416 getkyokaigap( egap1, seq1, cut1[i+1], clus1 );
1417 getkyokaigap( egap2, seq2, cut2[i+1], clus2 );
1421 for( j=0; j<clus1; j++ ) egap1[j] = 'o';
1422 for( j=0; j<clus2; j++ ) egap2[j] = 'o';
1426 fprintf( stderr, "kyokkaigap1(%d)=", cut1[i]-1 );
1427 for( j=0; j<clus1; j++ )
1428 fprintf( stderr, "%c", sgap1[j] );
1429 fprintf( stderr, "=kyokkaigap1-start\n" );
1432 fprintf( stderr, "kyokkaigap2(%d)=", cut2[i]-1 );
1433 for( j=0; j<clus2; j++ )
1434 fprintf( stderr, "%c", sgap2[j] );
1435 fprintf( stderr, "=kyokkaigap2-start\n" );
1438 fprintf( stderr, "kyokkaigap1(%d)=", cut1[i]-1 );
1439 for( j=0; j<clus1; j++ )
1440 fprintf( stderr, "%c", egap1[j] );
1441 fprintf( stderr, "=kyokkaigap1-end\n" );
1444 fprintf( stderr, "kyokkaigap2(%d)=", cut2[i]-1 );
1445 for( j=0; j<clus2; j++ )
1446 fprintf( stderr, "%c", egap2[j] );
1447 fprintf( stderr, "=kyokkaigap2-end\n" );
1452 fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
1455 fprintf( stderr, "DP %03d / %03d\r", i+1, count-1 );
1458 for( j=0; j<clus1; j++ )
1460 strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
1461 tmpres1[j][cut1[i+1]-cut1[i]] = 0;
1463 if( kobetsubunkatsu && fftkeika ) commongappick( clus1, tmpres1 ); //dvtditr
\e$B$K8F$P$l$?$H$-
\e(B fftkeika=1
1464 // if( kobetsubunkatsu ) commongappick( clus1, tmpres1 );
1465 for( j=0; j<clus2; j++ )
1467 strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
1468 tmpres2[j][cut2[i+1]-cut2[i]] = 0;
1470 if( kobetsubunkatsu && fftkeika ) commongappick( clus2, tmpres2 ); //dvtditr
\e$B$K8F$P$l$?$H$-
\e(B fftkeika=1
1471 // if( kobetsubunkatsu ) commongappick( clus2, tmpres2 );
1475 fprintf( stderr, "Not supported\n" );
1479 fprintf( stderr, "i=%d, before alignment", i );
1480 fprintf( stderr, "%4d\n", totallen );
1481 fprintf( stderr, "\n\n" );
1482 for( j=0; j<clus1; j++ )
1484 fprintf( stderr, "%s\n", tmpres1[j] );
1486 fprintf( stderr, "-------\n" );
1487 for( j=0; j<clus2; j++ )
1489 fprintf( stderr, "%s\n", tmpres2[j] );
1494 fprintf( stdout, "writing input\n" );
1495 for( j=0; j<clus1; j++ )
1497 fprintf( stdout, ">%d of GROUP1\n", j );
1498 fprintf( stdout, "%s\n", tmpres1[j] );
1500 for( j=0; j<clus2; j++ )
1502 fprintf( stdout, ">%d of GROUP2\n", j );
1503 fprintf( stdout, "%s\n", tmpres2[j] );
1510 totalscore += Aalign( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen );
1513 totalscore += MSalignmm( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres, headgp, tailgp );
1516 if( clus1 == 1 && clus2 == 1 )
1518 totalscore += G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
1521 totalscore += A__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres, headgp, tailgp );
1524 if( clus1 == 1 && clus2 == 1 )
1526 totalscore += G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
1529 totalscore += H__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, sgap1, sgap2, egap1, egap2 );
1532 if( clus1 == 1 && clus2 == 1 )
1534 totalscore += G__align11( tmpres1, tmpres2, alloclen, headgp, tailgp );
1537 totalscore += Q__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, NULL, &dumfl, sgap1, sgap2, egap1, egap2 );
1540 fprintf( stderr, "alg = %c\n", alg );
1541 ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
1545 #ifdef enablemultithread
1546 if( chudanres && *chudanres )
1548 // fprintf( stderr, "\n\n## CHUUDAN!!! at Falign_localhom\n" );
1553 nlen = strlen( tmpres1[0] );
1554 if( totallen + nlen > alloclen )
1556 fprintf( stderr, "totallen=%d + nlen=%d > alloclen = %d\n", totallen, nlen, alloclen );
1557 ErrorExit( "LENGTH OVER in Falign\n " );
1559 for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
1560 for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
1563 fprintf( stderr, "i=%d", i );
1564 fprintf( stderr, "%4d\n", totallen );
1565 fprintf( stderr, "\n\n" );
1566 for( j=0; j<clus1; j++ )
1568 fprintf( stderr, "%s\n", tmpres1[j] );
1570 fprintf( stderr, "-------\n" );
1571 for( j=0; j<clus2; j++ )
1573 fprintf( stderr, "%s\n", tmpres2[j] );
1578 fprintf( stderr, "DP ... done \n" );
1581 for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
1582 for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
1584 for( j=0; j<clus1; j++ )
1586 fprintf( stderr, "%s\n", result1[j] );
1588 fprintf( stderr, "- - - - - - - - - - -\n" );
1589 for( j=0; j<clus2; j++ )
1591 fprintf( stderr, "%s\n", result2[j] );
1594 return( totalscore );
1605 sakujo wo kentou (2010/10/05)
1607 float Falign_udpari_long( char **seq1, char **seq2,
1608 double *eff1, double *eff2,
1609 int clus1, int clus2,
1610 int alloclen, int *fftlog )
1612 int i, j, k, l, m, maxk;
1613 int nlen, nlen2, nlen4;
1614 static TLS int prevalloclen = 0;
1615 static TLS int crossscoresize = 0;
1616 static TLS char **tmpseq1 = NULL;
1617 static TLS char **tmpseq2 = NULL;
1618 static TLS char **tmpptr1 = NULL;
1619 static TLS char **tmpptr2 = NULL;
1620 static TLS char **tmpres1 = NULL;
1621 static TLS char **tmpres2 = NULL;
1622 static TLS char **result1 = NULL;
1623 static TLS char **result2 = NULL;
1625 static TLS char **rndseq1 = NULL;
1626 static TLS char **rndseq2 = NULL;
1628 static TLS Fukusosuu **seqVector1 = NULL;
1629 static TLS Fukusosuu **seqVector2 = NULL;
1630 static TLS Fukusosuu **naiseki = NULL;
1631 static TLS Fukusosuu *naisekiNoWa = NULL;
1632 static TLS double *soukan = NULL;
1633 static TLS double **crossscore = NULL;
1635 static TLS int *kouho = NULL;
1636 static TLS Segment *segment = NULL;
1637 static TLS Segment *segment1 = NULL;
1638 static TLS Segment *segment2 = NULL;
1639 static TLS Segment **sortedseg1 = NULL;
1640 static TLS Segment **sortedseg2 = NULL;
1641 static TLS int *cut1 = NULL;
1642 static TLS int *cut2 = NULL;
1643 static TLS char *sgap1, *egap1, *sgap2, *egap2;
1644 static TLS int localalloclen = 0;
1653 // float dumfl = 0.0;
1657 len1 = strlen( seq1[0] );
1658 len2 = strlen( seq2[0] );
1659 nlentmp = MAX( len1, len2 );
1662 while( nlentmp >= nlen ) nlen <<= 1;
1664 fprintf( stderr, "### nlen = %d\n", nlen );
1667 nlen2 = nlen/2; nlen4 = nlen2 / 2;
1670 fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
1671 fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
1674 if( prevalloclen != alloclen ) // Falign_noudp mo kaeru
1678 FreeCharMtx( result1 );
1679 FreeCharMtx( result2 );
1680 FreeCharMtx( tmpres1 );
1681 FreeCharMtx( tmpres2 );
1683 // fprintf( stderr, "\n\n\nreallocating ...\n" );
1684 result1 = AllocateCharMtx( njob, alloclen );
1685 result2 = AllocateCharMtx( njob, alloclen );
1686 tmpres1 = AllocateCharMtx( njob, alloclen );
1687 tmpres2 = AllocateCharMtx( njob, alloclen );
1688 prevalloclen = alloclen;
1691 if( !localalloclen )
1693 sgap1 = AllocateCharVec( njob );
1694 egap1 = AllocateCharVec( njob );
1695 sgap2 = AllocateCharVec( njob );
1696 egap2 = AllocateCharVec( njob );
1697 kouho = AllocateIntVec( NKOUHO_LONG );
1698 cut1 = AllocateIntVec( MAXSEG );
1699 cut2 = AllocateIntVec( MAXSEG );
1700 tmpptr1 = AllocateCharMtx( njob, 0 );
1701 tmpptr2 = AllocateCharMtx( njob, 0 );
1702 segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
1703 segment1 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
1704 segment2 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
1705 sortedseg1 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
1706 sortedseg2 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
1707 if( !( segment && segment1 && segment2 && sortedseg1 && sortedseg2 ) )
1708 ErrorExit( "Allocation error\n" );
1710 if ( scoremtx == -1 ) n20or4or2 = 1;
1711 else if( fftscore ) n20or4or2 = 1;
1712 else n20or4or2 = 20;
1714 if( localalloclen < nlen )
1719 if( !kobetsubunkatsu )
1721 FreeFukusosuuMtx ( seqVector1 );
1722 FreeFukusosuuMtx ( seqVector2 );
1723 FreeFukusosuuVec( naisekiNoWa );
1724 FreeFukusosuuMtx( naiseki );
1725 FreeDoubleVec( soukan );
1727 FreeCharMtx( tmpseq1 );
1728 FreeCharMtx( tmpseq2 );
1731 FreeCharMtx( rndseq1 );
1732 FreeCharMtx( rndseq2 );
1737 tmpseq1 = AllocateCharMtx( njob, nlen );
1738 tmpseq2 = AllocateCharMtx( njob, nlen );
1739 if( !kobetsubunkatsu )
1741 naisekiNoWa = AllocateFukusosuuVec( nlen );
1742 naiseki = AllocateFukusosuuMtx( n20or4or2, nlen );
1743 seqVector1 = AllocateFukusosuuMtx( n20or4or2, nlen+1 );
1744 seqVector2 = AllocateFukusosuuMtx( n20or4or2, nlen+1 );
1745 soukan = AllocateDoubleVec( nlen+1 );
1748 rndseq1 = AllocateCharMtx( njob, nlen );
1749 rndseq2 = AllocateCharMtx( njob, nlen );
1750 for( i=0; i<njob; i++ )
1752 generateRndSeq( rndseq1[i], nlen );
1753 generateRndSeq( rndseq2[i], nlen );
1756 localalloclen = nlen;
1759 for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
1760 for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );
1763 fftfp = fopen( "input_of_Falign", "w" );
1764 fprintf( fftfp, "nlen = %d\n", nlen );
1765 fprintf( fftfp, "seq1: ( %d sequences ) \n", clus1 );
1766 for( i=0; i<clus1; i++ )
1767 fprintf( fftfp, "%s\n", seq1[i] );
1768 fprintf( fftfp, "seq2: ( %d sequences ) \n", clus2 );
1769 for( i=0; i<clus2; i++ )
1770 fprintf( fftfp, "%s\n", seq2[i] );
1772 system( "less input_of_Falign < /dev/tty > /dev/tty" );
1774 if( !kobetsubunkatsu )
1776 fprintf( stderr, " FFT ... " );
1778 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
1779 if( scoremtx == -1 )
1781 for( i=0; i<clus1; i++ )
1782 seq_vec_4( seqVector1[0], eff1[i], tmpseq1[i] );
1786 for( i=0; i<clus1; i++ )
1789 seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
1790 seq_vec_2( seqVector1[1], volume, eff1[i], tmpseq1[i] );
1792 seq_vec_5( seqVector1[0], polarity, volume, eff1[i], tmpseq1[i] );
1798 for( i=0; i<clus1; i++ )
1799 seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
1802 for( i=0; i<clus1; i++ )
1804 vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
1808 fftfp = fopen( "seqVec", "w" );
1809 fprintf( fftfp, "before transform\n" );
1810 for( k=0; k<n20or4or2; k++ )
1812 fprintf( fftfp, "nlen=%d\n", nlen );
1813 fprintf( fftfp, "%c\n", amino[k] );
1814 for( l=0; l<nlen; l++ )
1815 fprintf( fftfp, "%f %f\n", seqVector1[k][l].R, seqVector1[k][l].I );
1818 system( "less seqVec < /dev/tty > /dev/tty" );
1821 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
1822 if( scoremtx == -1 )
1824 for( i=0; i<clus2; i++ )
1825 seq_vec_4( seqVector2[0], eff2[i], tmpseq2[i] );
1829 for( i=0; i<clus2; i++ )
1832 seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
1833 seq_vec_2( seqVector2[1], volume, eff2[i], tmpseq2[i] );
1835 seq_vec_5( seqVector2[0], polarity, volume, eff2[i], tmpseq2[i] );
1841 for( i=0; i<clus2; i++ )
1842 seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
1845 for( i=0; i<clus2; i++ )
1847 vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
1852 fftfp = fopen( "seqVec2", "w" );
1853 fprintf( fftfp, "before fft\n" );
1854 for( k=0; k<n20or4or2; k++ )
1856 fprintf( fftfp, "%c\n", amino[k] );
1857 for( l=0; l<nlen; l++ )
1858 fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
1861 system( "less seqVec2 < /dev/tty > /dev/tty" );
1864 for( j=0; j<n20or4or2; j++ )
1866 fft( nlen, seqVector2[j], 0 );
1867 fft( nlen, seqVector1[j], 0 );
1870 fftfp = fopen( "seqVec2", "w" );
1871 fprintf( fftfp, "#after fft\n" );
1872 for( k=0; k<n20or4or2; k++ )
1874 fprintf( fftfp, "#%c\n", amino[k] );
1875 for( l=0; l<nlen; l++ )
1876 fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
1879 system( "less seqVec2 < /dev/tty > /dev/tty" );
1882 for( k=0; k<n20or4or2; k++ )
1884 for( l=0; l<nlen; l++ )
1885 calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
1887 for( l=0; l<nlen; l++ )
1889 naisekiNoWa[l].R = 0.0;
1890 naisekiNoWa[l].I = 0.0;
1891 for( k=0; k<n20or4or2; k++ )
1893 naisekiNoWa[l].R += naiseki[k][l].R;
1894 naisekiNoWa[l].I += naiseki[k][l].I;
1899 fftfp = fopen( "naisekiNoWa", "w" );
1900 fprintf( fftfp, "#Before fft\n" );
1901 for( l=0; l<nlen; l++ )
1902 fprintf( fftfp, "%d %f %f\n", l, naisekiNoWa[l].R, naisekiNoWa[l].I );
1904 system( "less naisekiNoWa < /dev/tty > /dev/tty " );
1907 fft( -nlen, naisekiNoWa, 0 );
1909 for( m=0; m<=nlen2; m++ )
1910 soukan[m] = naisekiNoWa[nlen2-m].R;
1911 for( m=nlen2+1; m<nlen; m++ )
1912 soukan[m] = naisekiNoWa[nlen+nlen2-m].R;
1915 fftfp = fopen( "naisekiNoWa", "w" );
1916 fprintf( fftfp, "#After fft\n" );
1917 for( l=0; l<nlen; l++ )
1918 fprintf( fftfp, "%d %f\n", l, naisekiNoWa[l].R );
1920 fftfp = fopen( "list.plot", "w" );
1921 fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
1923 system( "/usr/bin/gnuplot list.plot &" );
1926 fprintf( stderr, "soukan\n" );
1927 for( l=0; l<nlen; l++ )
1928 fprintf( stderr, "%d %f\n", l-nlen2, soukan[l] );
1930 fftfp = fopen( "list.plot", "w" );
1931 fprintf( fftfp, "plot 'frt'\n pause +1" );
1933 system( "/usr/bin/gnuplot list.plot" );
1938 nkouho = getKouho( kouho, NKOUHO_LONG, soukan, nlen );
1941 for( i=0; i<nkouho; i++ )
1943 fprintf( stderr, "kouho[%d] = %d\n", i, kouho[i] );
1949 fprintf( stderr, "Searching anchors ... " );
1957 fftfp = fopen( "cand", "w" );
1960 if( kobetsubunkatsu )
1970 for( k=0; k<maxk; k++ )
1973 if( lag <= -len1 || len2 <= lag ) continue;
1974 // fprintf( stderr, "k=%d, lag=%d\n", k, lag );
1975 zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
1977 fftfp = fopen( "cand", "a" );
1978 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
1979 fprintf( fftfp, "%s\n", tmpptr1[0] );
1980 fprintf( fftfp, ">Candidate No.%d lag = %d\n", k+1, lag );
1981 fprintf( fftfp, "%s\n", tmpptr2[0] );
1982 fprintf( fftfp, ">\n", k+1, lag );
1986 // fprintf( stderr, "lag = %d\n", lag );
1987 tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );
1988 // fprintf( stderr, "lag = %d, %d found\n", lag, tmpint );
1990 // if( lag == -50 ) exit( 1 );
1992 if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );
1994 // fprintf( stderr, "##### k=%d / %d\n", k, maxk );
1995 // if( tmpint == 0 ) break; // 060430 iinoka ? // 090530 yameta
1996 while( tmpint-- > 0 )
1999 if( segment[count].end - segment[count].start < fftWinSize )
2007 segment1[count].start = segment[count].start ;
2008 segment1[count].end = segment[count].end ;
2009 segment1[count].center = segment[count].center;
2010 segment1[count].score = segment[count].score;
2012 segment2[count].start = segment[count].start + lag;
2013 segment2[count].end = segment[count].end + lag;
2014 segment2[count].center = segment[count].center + lag;
2015 segment2[count].score = segment[count].score ;
2019 segment1[count].start = segment[count].start - lag;
2020 segment1[count].end = segment[count].end - lag;
2021 segment1[count].center = segment[count].center - lag;
2022 segment1[count].score = segment[count].score ;
2024 segment2[count].start = segment[count].start ;
2025 segment2[count].end = segment[count].end ;
2026 segment2[count].center = segment[count].center;
2027 segment2[count].score = segment[count].score ;
2030 fprintf( stderr, "##### k=%d / %d\n", k, maxk );
2031 fprintf( stderr, "anchor %d, score = %f\n", count, segment1[count].score );
2032 fprintf( stderr, "in 1 %d\n", segment1[count].center );
2033 fprintf( stderr, "in 2 %d\n", segment2[count].center );
2035 segment1[count].pair = &segment2[count];
2036 segment2[count].pair = &segment1[count];
2039 fprintf( stderr, "count=%d\n", count );
2044 if( !kobetsubunkatsu )
2045 fprintf( stderr, "done. (%d anchors) ", count );
2047 if( !count && fftNoAnchStop )
2048 ErrorExit( "Cannot detect anchor!" );
2050 fprintf( stderr, "RESULT before sort:\n" );
2051 for( l=0; l<count+1; l++ )
2053 fprintf( stderr, "cut[%d]=%d, ", l, segment1[l].center );
2054 fprintf( stderr, "%d score = %f\n", segment2[l].center, segment1[l].score );
2058 for( i=0; i<count; i++ )
2060 sortedseg1[i] = &segment1[i];
2061 sortedseg2[i] = &segment2[i];
2064 tmpsort( count, sortedseg1 );
2065 tmpsort( count, sortedseg2 );
2066 qsort( sortedseg1, count, sizeof( Segment * ), segcmp );
2067 qsort( sortedseg2, count, sizeof( Segment * ), segcmp );
2069 mymergesort( 0, count-1, sortedseg1 );
2070 mymergesort( 0, count-1, sortedseg2 );
2072 for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
2073 for( i=0; i<count; i++ ) sortedseg2[i]->number = i;
2077 if( kobetsubunkatsu )
2079 for( i=0; i<count; i++ )
2081 cut1[i+1] = sortedseg1[i]->center;
2082 cut2[i+1] = sortedseg2[i]->center;
2086 cut1[count+1] = len1;
2087 cut2[count+1] = len2;
2095 if( crossscoresize < count+2 )
2097 crossscoresize = count+2;
2099 if( fftkeika ) fprintf( stderr, "######allocating crossscore, size = %d\n", crossscoresize );
2101 if( crossscore ) FreeDoubleMtx( crossscore );
2102 crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
2104 for( i=0; i<count+2; i++ ) for( j=0; j<count+2; j++ )
2105 crossscore[i][j] = 0.0;
2106 for( i=0; i<count; i++ )
2108 crossscore[segment1[i].number+1][segment1[i].pair->number+1] = segment1[i].score;
2109 cut1[i+1] = sortedseg1[i]->center;
2110 cut2[i+1] = sortedseg2[i]->center;
2114 fprintf( stderr, "AFTER SORT\n" );
2115 for( i=0; i<count+1; i++ ) fprintf( stderr, "%d, %d\n", cut1[i], cut2[i] );
2116 fprintf( stderr, "crossscore = \n" );
2117 for( i=0; i<count+1; i++ )
2119 for( j=0; j<count+1; j++ )
2120 fprintf( stderr, "%.0f ", crossscore[i][j] );
2121 fprintf( stderr, "\n" );
2125 crossscore[0][0] = 10000000.0;
2128 crossscore[count+1][count+1] = 10000000.0;
2129 cut1[count+1] = len1;
2130 cut2[count+1] = len2;
2134 // fprintf( stderr, "\n\n\ncalling blockAlign2\n\n\n\n" );
2135 blockAlign2( cut1, cut2, sortedseg1, sortedseg2, crossscore, &count );
2137 // if( count-count0 )
2138 // fprintf( stderr, "%d unused anchors\n", count0-count );
2140 if( !kobetsubunkatsu && fftkeika )
2141 fprintf( stderr, "%d anchors found\n", count );
2144 if( count0 > count )
2147 fprintf( stderr, "\7 REPEAT!? \n" );
2149 fprintf( stderr, "REPEAT!? \n" );
2151 if( fftRepeatStop ) exit( 1 );
2154 else fprintf( stderr, "done\n" );
2162 fprintf( stderr, "\nMany anchors were found. The upper-level DP is skipped.\n\n" );
2167 for( i=0; i<count; i++ )
2169 // fprintf( stderr, "i=%d, %d-%d ?\n", i, sortedseg1[i]->center, sortedseg1[i]->pair->center );
2170 if( sortedseg1[i]->center > cut1[count0]
2171 && sortedseg1[i]->pair->center > cut2[count0] )
2174 cut1[count0] = sortedseg1[i]->center;
2175 cut2[count0] = sortedseg1[i]->pair->center;
2179 if( i && sortedseg1[i]->score > sortedseg1[i-1]->score )
2181 if( sortedseg1[i]->center > cut1[count0-1]
2182 && sortedseg1[i]->pair->center > cut2[count0-1] )
2184 cut1[count0] = sortedseg1[i]->center;
2185 cut2[count0] = sortedseg1[i]->pair->center;
2194 // if( count-count0 )
2195 // fprintf( stderr, "%d anchors unused\n", count-count0 );
2196 cut1[count0+1] = len1;
2197 cut2[count0+1] = len2;
2207 fftfp = fopen( "fft", "a" );
2208 fprintf( fftfp, "RESULT after sort:\n" );
2209 for( l=0; l<count; l++ )
2211 fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
2212 fprintf( fftfp, "%d\n", segment2[l].center );
2218 fprintf( stderr, "RESULT after blckalign:\n" );
2219 for( l=0; l<count+1; l++ )
2221 fprintf( stderr, "cut : %d %d\n", cut1[l], cut2[l] );
2226 fprintf( trap_g, "Devided to %d segments\n", count-1 );
2227 fprintf( trap_g, "%d %d forg\n", MIN( clus1, clus2 ), count-1 );
2231 for( j=0; j<clus1; j++ ) result1[j][0] = 0;
2232 for( j=0; j<clus2; j++ ) result2[j][0] = 0;
2235 for( i=0; i<count-1; i++ )
2238 if( i == 0 ) headgp = outgap; else headgp = 1;
2239 if( i == count-2 ) tailgp = outgap; else tailgp = 1;
2243 // getkyokaigap( sgap1, seq1, cut1[i]-1, clus1 );
2244 // getkyokaigap( sgap2, seq2, cut2[i]-1, clus2 );
2245 getkyokaigap( sgap1, tmpres1, nlen-1, clus1 );
2246 getkyokaigap( sgap2, tmpres2, nlen-1, clus2 );
2250 for( j=0; j<clus1; j++ ) sgap1[j] = 'o';
2251 for( j=0; j<clus2; j++ ) sgap2[j] = 'o';
2253 if( cut1[i+1] != len1 )
2255 getkyokaigap( egap1, seq1, cut1[i+1], clus1 );
2256 getkyokaigap( egap2, seq2, cut2[i+1], clus2 );
2260 for( j=0; j<clus1; j++ ) egap1[j] = 'o';
2261 for( j=0; j<clus2; j++ ) egap2[j] = 'o';
2264 fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
2267 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 );
2270 for( j=0; j<clus1; j++ )
2272 strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
2273 tmpres1[j][cut1[i+1]-cut1[i]] = 0;
2275 if( kobetsubunkatsu && fftkeika ) commongappick( clus1, tmpres1 ); //dvtditr
\e$B$K8F$P$l$?$H$-
\e(B fftkeika=1
2276 // if( kobetsubunkatsu ) commongappick( clus1, tmpres1 );
2277 for( j=0; j<clus2; j++ )
2279 // fprintf( stderr, "### cut2[i+1]-cut2[i] = %d\n", cut2[i+1]-cut2[i] );
2280 if( cut2[i+1]-cut2[i] <= 0 )
2281 fprintf( stderr, "### cut2[i+1]=%d, cut2[i]=%d\n", cut2[i+1], cut2[i] );
2282 strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
2283 tmpres2[j][cut2[i+1]-cut2[i]] = 0;
2285 if( kobetsubunkatsu && fftkeika ) commongappick( clus2, tmpres2 ); //dvtditr
\e$B$K8F$P$l$?$H$-
\e(B fftkeika=1
2286 // if( kobetsubunkatsu ) commongappick( clus2, tmpres2 );
2290 fprintf( stderr, "Not supported\n" );
2294 fprintf( stderr, "i=%d, before alignment", i );
2295 fprintf( stderr, "%4d\n", totallen );
2296 fprintf( stderr, "\n\n" );
2297 for( j=0; j<clus1; j++ )
2299 fprintf( stderr, "%s\n", tmpres1[j] );
2301 fprintf( stderr, "-------\n" );
2302 for( j=0; j<clus2; j++ )
2304 fprintf( stderr, "%s\n", tmpres2[j] );
2309 fprintf( stdout, "writing input\n" );
2310 for( j=0; j<clus1; j++ )
2312 fprintf( stdout, ">%d of GROUP1\n", j );
2313 fprintf( stdout, "%s\n", tmpres1[j] );
2315 for( j=0; j<clus2; j++ )
2317 fprintf( stdout, ">%d of GROUP2\n", j );
2318 fprintf( stdout, "%s\n", tmpres2[j] );
2325 totalscore += MSalignmm( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, sgap1, sgap2, egap1, egap2, NULL, 0, NULL, headgp, tailgp );
2328 fprintf( stderr, "alg = %c\n", alg );
2329 ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
2333 nlen = strlen( tmpres1[0] );
2334 if( totallen + nlen > alloclen )
2336 fprintf( stderr, "totallen=%d + nlen=%d > alloclen = %d\n", totallen, nlen, alloclen );
2337 ErrorExit( "LENGTH OVER in Falign\n " );
2339 for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
2340 for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
2343 fprintf( stderr, "i=%d", i );
2344 fprintf( stderr, "%4d\n", totallen );
2345 fprintf( stderr, "\n\n" );
2346 for( j=0; j<clus1; j++ )
2348 fprintf( stderr, "%s\n", tmpres1[j] );
2350 fprintf( stderr, "-------\n" );
2351 for( j=0; j<clus2; j++ )
2353 fprintf( stderr, "%s\n", tmpres2[j] );
2358 fprintf( stderr, "DP ... done \n" );
2361 for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
2362 for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
2364 for( j=0; j<clus1; j++ )
2366 fprintf( stderr, "%s\n", result1[j] );
2368 fprintf( stderr, "- - - - - - - - - - -\n" );
2369 for( j=0; j<clus2; j++ )
2371 fprintf( stderr, "%s\n", result2[j] );
2374 return( totalscore );