1 #include "mltaln.h"
2
3 //static FILE *fftfp;
4 static TLS int n20or4or2;
5
6 #define KEIKA 0
7 #define RND 0
8 #define DEBUG 0
9
10 extern int fft( int, Fukusosuu *, int );
11
12
13 #if 0
14 static void generateRndSeq( char *seq, int len )
15 {
16 while( len-- )
17 #if 1
18 *seq++ = (int)( rnd() * n20or4or2 );
19 #else
20 *seq++ = (int)1;
21 #endif
22 }
23 #endif
24
vec_init(Fukusosuu * result,int nlen)25 static void vec_init( Fukusosuu *result, int nlen )
26 {
27 while( nlen-- )
28 {
29 result->R = result->I = 0.0;
30 result++;
31 }
32 }
33
34 #if 0
35 static void vec_init2( Fukusosuu **result, char *seq, double eff, int st, int ed )
36 {
37 int i;
38 for( i=st; i<ed; i++ )
39 result[(int)*seq++][i].R += eff;
40 }
41 #endif
42
seq_vec_2(Fukusosuu * result,double * score,double incr,char * seq)43 static void seq_vec_2( Fukusosuu *result, double *score, double incr, char *seq )
44 {
45 static TLS int n;
46 for( ; *seq; result++ )
47 {
48 n = amino_n[(int)*seq++];
49 if( n < 20 && n >= 0 ) result->R += incr * score[n];
50 #if 0
51 fprintf( stderr, "n=%d, score=%f, inc=%f R=%f\n",n, score[n], incr * score[n], result->R );
52 #endif
53 }
54 }
55
seq_vec_3(Fukusosuu ** result,double incr,char * seq)56 static void seq_vec_3( Fukusosuu **result, double incr, char *seq )
57 {
58 int i;
59 int n;
60 for( i=0; *seq; i++ )
61 {
62 n = amino_n[(int)*seq++];
63 if( n < n20or4or2 && n >= 0 ) result[n][i].R += incr;
64 }
65 }
66
67
68 #if 0
69 static void seq_vec( Fukusosuu *result, char query, double incr, char *seq )
70 {
71 #if 0
72 int bk = nlen;
73 #endif
74 while( *seq )
75 {
76 if( *seq++ == query ) result->R += incr;
77 result++;
78 #if 0
79 fprintf( stderr, "i = %d result->R = %f\n", bk-nlen, (result-1)->R );
80 #endif
81 }
82 }
83
84 static int checkRepeat( int num, int *cutpos )
85 {
86 int tmp, buf;
87
88 buf = *cutpos;
89 while( num-- )
90 {
91 if( ( tmp = *cutpos++ ) < buf ) return( 1 );
92 buf = tmp;
93 }
94 return( 0 );
95 }
96
97 static int segcmp( void *ptr1, void *ptr2 )
98 {
99 int diff;
100 Segment **seg1 = (Segment **)ptr1;
101 Segment **seg2 = (Segment **)ptr2;
102 #if 0
103 return( (*seg1)->center - (*seg2)->center );
104 #else
105 diff = (*seg1)->center - (*seg2)->center;
106 if( diff ) return( diff );
107
108 diff = (*seg1)->start - (*seg2)->start;
109 if( diff ) return( diff );
110
111 diff = (*seg1)->end - (*seg2)->end;
112 if( diff ) return( diff );
113
114 fprintf( stderr, "USE STABLE SORT !!\n" );
115 exit( 1 );
116 return( 0 );
117 #endif
118 }
119
120 #endif
121
122
mymergesort(int first,int last,Segment ** seg)123 static void mymergesort( int first, int last, Segment **seg )
124 {
125 int middle;
126 static TLS int i, j, k, p;
127 static TLS int allo = 0;
128 static TLS Segment **work = NULL;
129
130 if( seg == NULL )
131 {
132 free( work ); work = NULL;
133 return;
134 }
135
136 if( last > allo )
137 {
138 allo = last;
139 if( work ) free( work );
140 work = (Segment **)calloc( allo / 2 + 1, sizeof( Segment *) );
141 }
142
143 if( first < last )
144 {
145 middle = ( first + last ) / 2;
146 mymergesort( first, middle, seg );
147 mymergesort( middle+1, last, seg );
148 p = 0;
149 for( i=first; i<=middle; i++ ) work[p++] = seg[i];
150 i = middle + 1; j = 0; k = first;
151 while( i <= last && j < p )
152 {
153 if( work[j]->center <= seg[i]->center )
154 seg[k++] = work[j++];
155 else
156 seg[k++] = seg[i++];
157 }
158 while( j < p ) seg[k++] = work[j++];
159 }
160 }
161
162
Falign_localhom(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,LocalHom *** localhom,double * totalimpmatch,int * gapmap1,int * gapmap2,int * chudanpt,int chudanref,int * chudanres)163 double Falign_localhom( int **whichmtx, double ***scoringmatrices, double **n_dynamicmtx,
164 char **seq1, char **seq2,
165 double *eff1, double *eff2,
166 double **eff1s, double **eff2s,
167 int clus1, int clus2,
168 int alloclen,
169 LocalHom ***localhom, double *totalimpmatch,
170 int *gapmap1, int *gapmap2,
171 int *chudanpt, int chudanref, int *chudanres )
172 {
173 // tditeration.c deha alloclen ha huhen nanode
174 // prevalloclen ha iranai.
175 int i, j, k, l, m, maxk;
176 int nlen, nlen2, nlen4;
177 static TLS int crossscoresize = 0;
178 static TLS char **tmpseq1 = NULL;
179 static TLS char **tmpseq2 = NULL;
180 static TLS char **tmpptr1 = NULL;
181 static TLS char **tmpptr2 = NULL;
182 static TLS char **tmpres1 = NULL;
183 static TLS char **tmpres2 = NULL;
184 static TLS char **result1 = NULL;
185 static TLS char **result2 = NULL;
186 #if RND
187 static TLS char **rndseq1 = NULL;
188 static TLS char **rndseq2 = NULL;
189 #endif
190 static TLS Fukusosuu **seqVector1 = NULL;
191 static TLS Fukusosuu **seqVector2 = NULL;
192 static TLS Fukusosuu **naiseki = NULL;
193 static TLS Fukusosuu *naisekiNoWa = NULL;
194 static TLS double *soukan = NULL;
195 static TLS double **crossscore = NULL;
196 int nlentmp;
197 static TLS int *kouho = NULL;
198 static TLS Segment *segment = NULL;
199 static TLS Segment *segment1 = NULL;
200 static TLS Segment *segment2 = NULL;
201 static TLS Segment **sortedseg1 = NULL;
202 static TLS Segment **sortedseg2 = NULL;
203 static TLS int *cut1 = NULL;
204 static TLS int *cut2 = NULL;
205 static TLS char *sgap1, *egap1, *sgap2, *egap2;
206 static TLS int localalloclen = 0;
207 int lag;
208 int tmpint;
209 int count, count0;
210 int len1, len2;
211 int totallen;
212 double totalscore;
213 double impmatch;
214
215 extern Fukusosuu *AllocateFukusosuuVec();
216 extern Fukusosuu **AllocateFukusosuuMtx();
217
218 if( seq1 == NULL )
219 {
220 if( result1 )
221 {
222 // fprintf( stderr, "Freeing localarrays in Falign\n" );
223 localalloclen = 0;
224 crossscoresize = 0;
225 mymergesort( 0, 0, NULL );
226 alignableReagion( 0, 0, NULL, NULL, NULL, NULL, NULL );
227 fft( 0, NULL, 1 );
228 // A__align( NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 ); // iru?
229 G__align11( NULL, NULL, NULL, 0, 0, 0 );
230 partA__align( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL );
231 partA__align_variousdist( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL );
232 blockAlign2( NULL, NULL, NULL, NULL, NULL, NULL );
233 if( crossscore ) FreeDoubleMtx( crossscore );
234 FreeCharMtx( result1 );
235 FreeCharMtx( result2 );
236 FreeCharMtx( tmpres1 );
237 FreeCharMtx( tmpres2 );
238 FreeCharMtx( tmpseq1 );
239 FreeCharMtx( tmpseq2 );
240 free( sgap1 );
241 free( egap1 );
242 free( sgap2 );
243 free( egap2 );
244 free( kouho );
245 free( cut1 );
246 free( cut2 );
247 free( tmpptr1 );
248 free( tmpptr2 );
249 free( segment );
250 free( segment1 );
251 free( segment2 );
252 free( sortedseg1 );
253 free( sortedseg2 );
254 if( !kobetsubunkatsu )
255 {
256 FreeFukusosuuMtx ( seqVector1 );
257 FreeFukusosuuMtx ( seqVector2 );
258 FreeFukusosuuVec( naisekiNoWa );
259 FreeFukusosuuMtx( naiseki );
260 FreeDoubleVec( soukan );
261 }
262 }
263 else
264 {
265 // fprintf( stderr, "Did not allocate localarrays in Falign\n" );
266 }
267
268 return( 0.0 );
269 }
270
271 len1 = strlen( seq1[0] );
272 len2 = strlen( seq2[0] );
273 nlentmp = MAX( len1, len2 );
274
275 nlen = 1;
276 while( nlentmp >= nlen ) nlen <<= 1;
277 #if 0
278 fprintf( stderr, "### nlen = %d\n", nlen );
279 #endif
280
281 nlen2 = nlen/2; nlen4 = nlen2 / 2;
282
283 #if DEBUG
284 fprintf( stderr, "len1 = %d, len2 = %d\n", len1, len2 );
285 fprintf( stderr, "nlentmp = %d, nlen = %d\n", nlentmp, nlen );
286 #endif
287
288 if( !localalloclen )
289 {
290 sgap1 = AllocateCharVec( njob );
291 egap1 = AllocateCharVec( njob );
292 sgap2 = AllocateCharVec( njob );
293 egap2 = AllocateCharVec( njob );
294 kouho = AllocateIntVec( NKOUHO );
295 cut1 = AllocateIntVec( MAXSEG );
296 cut2 = AllocateIntVec( MAXSEG );
297 tmpptr1 = AllocateCharMtx( njob, 0 );
298 tmpptr2 = AllocateCharMtx( njob, 0 );
299 result1 = AllocateCharMtx( njob, alloclen );
300 result2 = AllocateCharMtx( njob, alloclen );
301 tmpres1 = AllocateCharMtx( njob, alloclen );
302 tmpres2 = AllocateCharMtx( njob, alloclen );
303 // crossscore = AllocateDoubleMtx( MAXSEG, MAXSEG );
304 segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
305 segment1 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
306 segment2 = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
307 sortedseg1 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
308 sortedseg2 = (Segment **)calloc( MAXSEG, sizeof( Segment * ) );
309 if( !( segment && segment1 && segment2 && sortedseg1 && sortedseg2 ) )
310 ErrorExit( "Allocation error\n" );
311
312 if ( scoremtx == -1 ) n20or4or2 = 4;
313 else if( fftscore == 1 ) n20or4or2 = 2;
314 else n20or4or2 = 20;
315 }
316 if( localalloclen < nlen )
317 {
318 if( localalloclen )
319 {
320 #if 1
321 if( !kobetsubunkatsu )
322 {
323 FreeFukusosuuMtx ( seqVector1 );
324 FreeFukusosuuMtx ( seqVector2 );
325 FreeFukusosuuVec( naisekiNoWa );
326 FreeFukusosuuMtx( naiseki );
327 FreeDoubleVec( soukan );
328 }
329 FreeCharMtx( tmpseq1 );
330 FreeCharMtx( tmpseq2 );
331 #endif
332 #if RND
333 FreeCharMtx( rndseq1 );
334 FreeCharMtx( rndseq2 );
335 #endif
336 }
337
338 tmpseq1 = AllocateCharMtx( njob, nlen );
339 tmpseq2 = AllocateCharMtx( njob, nlen );
340 if( !kobetsubunkatsu )
341 {
342 naisekiNoWa = AllocateFukusosuuVec( nlen );
343 naiseki = AllocateFukusosuuMtx( n20or4or2, nlen );
344 seqVector1 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
345 seqVector2 = AllocateFukusosuuMtx( n20or4or2+1, nlen+1 );
346 soukan = AllocateDoubleVec( nlen+1 );
347 }
348 #if RND
349 rndseq1 = AllocateCharMtx( njob, nlen );
350 rndseq2 = AllocateCharMtx( njob, nlen );
351 for( i=0; i<njob; i++ )
352 {
353 generateRndSeq( rndseq1[i], nlen );
354 generateRndSeq( rndseq2[i], nlen );
355 }
356 #endif
357 localalloclen = nlen;
358 }
359
360 for( j=0; j<clus1; j++ ) strcpy( tmpseq1[j], seq1[j] );
361 for( j=0; j<clus2; j++ ) strcpy( tmpseq2[j], seq2[j] );
362
363 #if 0
364 fftfp = fopen( "input_of_Falign", "w" );
365 fprintf( fftfp, "nlen = %d\n", nlen );
366 fprintf( fftfp, "seq1: ( %d sequences ) \n", clus1 );
367 for( i=0; i<clus1; i++ )
368 fprintf( fftfp, "%s\n", seq1[i] );
369 fprintf( fftfp, "seq2: ( %d sequences ) \n", clus2 );
370 for( i=0; i<clus2; i++ )
371 fprintf( fftfp, "%s\n", seq2[i] );
372 fclose( fftfp );
373 system( "less input_of_Falign < /dev/tty > /dev/tty" );
374 #endif
375 if( !kobetsubunkatsu )
376 {
377 fprintf( stderr, "FFT ... " );
378
379 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector1[j], nlen );
380 if( fftscore && scoremtx != -1 )
381 {
382 for( i=0; i<clus1; i++ )
383 {
384 seq_vec_2( seqVector1[0], polarity, eff1[i], tmpseq1[i] );
385 seq_vec_2( seqVector1[1], volume, eff1[i], tmpseq1[i] );
386 }
387 }
388 else
389 {
390 #if 0
391 for( i=0; i<clus1; i++ ) for( j=0; j<n20or4or2; j++ )
392 seq_vec( seqVector1[j], amino[j], eff1[i], tmpseq1[i] );
393 #else
394 for( i=0; i<clus1; i++ )
395 seq_vec_3( seqVector1, eff1[i], tmpseq1[i] );
396 #endif
397 }
398 #if RND
399 for( i=0; i<clus1; i++ )
400 {
401 vec_init2( seqVector1, rndseq1[i], eff1[i], len1, nlen );
402 }
403 #endif
404 #if 0
405 fftfp = fopen( "seqVec", "w" );
406 fprintf( fftfp, "before transform\n" );
407 for( k=0; k<n20or4or2; k++ )
408 {
409 fprintf( fftfp, "nlen=%d\n", nlen );
410 fprintf( fftfp, "%c\n", amino[k] );
411 for( l=0; l<nlen; l++ )
412 fprintf( fftfp, "%f %f\n", seqVector1[k][l].R, seqVector1[k][l].I );
413 }
414 fclose( fftfp );
415 system( "less seqVec < /dev/tty > /dev/tty" );
416 #endif
417
418 for( j=0; j<n20or4or2; j++ ) vec_init( seqVector2[j], nlen );
419 if( fftscore && scoremtx != -1 )
420 {
421 for( i=0; i<clus2; i++ )
422 {
423 seq_vec_2( seqVector2[0], polarity, eff2[i], tmpseq2[i] );
424 seq_vec_2( seqVector2[1], volume, eff2[i], tmpseq2[i] );
425 }
426 }
427 else
428 {
429 #if 0
430 for( i=0; i<clus2; i++ ) for( j=0; j<n20or4or2; j++ )
431 seq_vec( seqVector2[j], amino[j], eff2[i], tmpseq2[i] );
432 #else
433 for( i=0; i<clus2; i++ )
434 seq_vec_3( seqVector2, eff2[i], tmpseq2[i] );
435 #endif
436 }
437 #if RND
438 for( i=0; i<clus2; i++ )
439 {
440 vec_init2( seqVector2, rndseq2[i], eff2[i], len2, nlen );
441 }
442 #endif
443
444 #if 0
445 fftfp = fopen( "seqVec2", "w" );
446 fprintf( fftfp, "before fft\n" );
447 for( k=0; k<n20or4or2; k++ )
448 {
449 fprintf( fftfp, "%c\n", amino[k] );
450 for( l=0; l<nlen; l++ )
451 fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
452 }
453 fclose( fftfp );
454 system( "less seqVec2 < /dev/tty > /dev/tty" );
455 #endif
456
457 for( j=0; j<n20or4or2; j++ )
458 {
459 fft( nlen, seqVector2[j], (j==0) );
460 fft( nlen, seqVector1[j], 0 );
461 }
462 #if 0
463 fftfp = fopen( "seqVec2", "w" );
464 fprintf( fftfp, "#after fft\n" );
465 for( k=0; k<n20or4or2; k++ )
466 {
467 fprintf( fftfp, "#%c\n", amino[k] );
468 for( l=0; l<nlen; l++ )
469 fprintf( fftfp, "%f %f\n", seqVector2[k][l].R, seqVector2[k][l].I );
470 }
471 fclose( fftfp );
472 system( "less seqVec2 < /dev/tty > /dev/tty" );
473 #endif
474
475 for( k=0; k<n20or4or2; k++ )
476 {
477 for( l=0; l<nlen; l++ )
478 calcNaiseki( naiseki[k]+l, seqVector1[k]+l, seqVector2[k]+l );
479 }
480 for( l=0; l<nlen; l++ )
481 {
482 naisekiNoWa[l].R = 0.0;
483 naisekiNoWa[l].I = 0.0;
484 for( k=0; k<n20or4or2; k++ )
485 {
486 naisekiNoWa[l].R += naiseki[k][l].R;
487 naisekiNoWa[l].I += naiseki[k][l].I;
488 }
489 }
490
491 #if 0
492 fftfp = fopen( "naisekiNoWa", "w" );
493 fprintf( fftfp, "#Before fft\n" );
494 for( l=0; l<nlen; l++ )
495 fprintf( fftfp, "%d %f %f\n", l, naisekiNoWa[l].R, naisekiNoWa[l].I );
496 fclose( fftfp );
497 system( "less naisekiNoWa < /dev/tty > /dev/tty " );
498 #endif
499
500 fft( -nlen, naisekiNoWa, 0 );
501
502 for( m=0; m<=nlen2; m++ )
503 soukan[m] = naisekiNoWa[nlen2-m].R;
504 for( m=nlen2+1; m<nlen; m++ )
505 soukan[m] = naisekiNoWa[nlen+nlen2-m].R;
506
507 #if 0
508 fftfp = fopen( "naisekiNoWa", "w" );
509 fprintf( fftfp, "#After fft\n" );
510 for( l=0; l<nlen; l++ )
511 fprintf( fftfp, "%d %f\n", l, naisekiNoWa[l].R );
512 fclose( fftfp );
513 fftfp = fopen( "list.plot", "w" );
514 fprintf( fftfp, "plot 'naisekiNoWa'\npause -1" );
515 fclose( fftfp );
516 system( "/usr/bin/gnuplot list.plot &" );
517 #endif
518 #if 0
519 fprintf( stderr, "frt write start\n" );
520 fftfp = fopen( "frt", "w" );
521 for( l=0; l<nlen; l++ )
522 fprintf( fftfp, "%d %f\n", l-nlen2, soukan[l] );
523 fclose( fftfp );
524 system( "less frt < /dev/tty > /dev/tty" );
525 #if 0
526 fftfp = fopen( "list.plot", "w" );
527 fprintf( fftfp, "plot 'frt'\n pause +1" );
528 fclose( fftfp );
529 system( "/usr/bin/gnuplot list.plot" );
530 #endif
531 #endif
532
533
534 getKouho( kouho, NKOUHO, soukan, nlen );
535
536 #if 0
537 for( i=0; i<NKOUHO; i++ )
538 {
539 fprintf( stderr, "kouho[%d] = %d\n", i, kouho[i] );
540 }
541 #endif
542 }
543
544 #if KEIKA
545 fprintf( stderr, "Searching anchors ... " );
546 #endif
547 count = 0;
548
549
550
551 #define CAND 0
552 #if CAND
553 fftfp = fopen( "cand", "w" );
554 fclose( fftfp );
555 #endif
556 if( kobetsubunkatsu )
557 {
558 maxk = 1;
559 kouho[0] = 0;
560 }
561 else
562 {
563 maxk = NKOUHO;
564 }
565
566 for( k=0; k<maxk; k++ )
567 {
568 lag = kouho[k];
569 zurasu2( lag, clus1, clus2, seq1, seq2, tmpptr1, tmpptr2 );
570 #if CAND
571 fftfp = fopen( "cand", "a" );
572 fprintf( fftfp, "Candidate No.%d lag = %d\n", k+1, lag );
573 fprintf( fftfp, "%s\n", tmpptr1[0] );
574 fprintf( fftfp, "%s\n", tmpptr2[0] );
575 fclose( fftfp );
576 #endif
577 tmpint = alignableReagion( clus1, clus2, tmpptr1, tmpptr2, eff1, eff2, segment+count );
578
579 if( count+tmpint > MAXSEG -3 ) ErrorExit( "TOO MANY SEGMENTS.\n" );
580
581
582 while( tmpint-- > 0 )
583 {
584 if( lag > 0 )
585 {
586 segment1[count].start = segment[count].start ;
587 segment1[count].end = segment[count].end ;
588 segment1[count].center = segment[count].center;
589 segment1[count].score = segment[count].score;
590
591 segment2[count].start = segment[count].start + lag;
592 segment2[count].end = segment[count].end + lag;
593 segment2[count].center = segment[count].center + lag;
594 segment2[count].score = segment[count].score ;
595 }
596 else
597 {
598 segment1[count].start = segment[count].start - lag;
599 segment1[count].end = segment[count].end - lag;
600 segment1[count].center = segment[count].center - lag;
601 segment1[count].score = segment[count].score ;
602
603 segment2[count].start = segment[count].start ;
604 segment2[count].end = segment[count].end ;
605 segment2[count].center = segment[count].center;
606 segment2[count].score = segment[count].score ;
607 }
608 #if 0
609 fftfp = fopen( "cand", "a" );
610 fprintf( fftfp, "Goukaku=%dko\n", tmpint );
611 fprintf( fftfp, "in 1 %d\n", segment1[count].center );
612 fprintf( fftfp, "in 2 %d\n", segment2[count].center );
613 fclose( fftfp );
614 #endif
615 segment1[count].pair = &segment2[count];
616 segment2[count].pair = &segment1[count];
617 count++;
618 #if 0
619 fprintf( stderr, "count=%d\n", count );
620 #endif
621 }
622 }
623 #if 1
624 if( !kobetsubunkatsu )
625 fprintf( stderr, "%d segments found\n", count );
626 #endif
627 if( !count && fftNoAnchStop )
628 ErrorExit( "Cannot detect anchor!" );
629 #if 0
630 fftfp = fopen( "fft", "a" );
631 fprintf( fftfp, "RESULT before sort:\n" );
632 for( l=0; l<count; l++ )
633 {
634 fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
635 fprintf( fftfp, "%d score = %f\n", segment2[l].center, segment1[l].score );
636 }
637 fclose( fftfp );
638 #endif
639
640 #if KEIKA
641 fprintf( stderr, "Aligning anchors ... " );
642 #endif
643 for( i=0; i<count; i++ )
644 {
645 sortedseg1[i] = &segment1[i];
646 sortedseg2[i] = &segment2[i];
647 }
648 #if 0
649 tmpsort( count, sortedseg1 );
650 tmpsort( count, sortedseg2 );
651 qsort( sortedseg1, count, sizeof( Segment * ), segcmp );
652 qsort( sortedseg2, count, sizeof( Segment * ), segcmp );
653 #else
654 mymergesort( 0, count-1, sortedseg1 );
655 mymergesort( 0, count-1, sortedseg2 );
656 #endif
657 for( i=0; i<count; i++ ) sortedseg1[i]->number = i;
658 for( i=0; i<count; i++ ) sortedseg2[i]->number = i;
659
660
661 if( kobetsubunkatsu )
662 {
663 for( i=0; i<count; i++ )
664 {
665 cut1[i+1] = sortedseg1[i]->center;
666 cut2[i+1] = sortedseg2[i]->center;
667 }
668 cut1[0] = 0;
669 cut2[0] = 0;
670 cut1[count+1] = len1;
671 cut2[count+1] = len2;
672 count += 2;
673 }
674 else
675 {
676 if( crossscoresize < count+2 )
677 {
678 crossscoresize = count+2;
679 #if 1
680 fprintf( stderr, "######allocating crossscore, size = %d\n", crossscoresize );
681 #endif
682 if( crossscore ) FreeDoubleMtx( crossscore );
683 crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
684 }
685 for( i=0; i<count+2; i++ ) for( j=0; j<count+2; j++ )
686 crossscore[i][j] = 0.0;
687 for( i=0; i<count; i++ )
688 {
689 crossscore[segment1[i].number+1][segment1[i].pair->number+1] = segment1[i].score;
690 cut1[i+1] = sortedseg1[i]->center;
691 cut2[i+1] = sortedseg2[i]->center;
692 }
693
694 #if DEBUG
695 fprintf( stderr, "AFTER SORT\n" );
696 for( i=0; i<count; i++ ) fprintf( stderr, "%d, %d\n", segment1[i].start, segment2[i].start );
697 #endif
698
699 crossscore[0][0] = 10000000.0;
700 cut1[0] = 0;
701 cut2[0] = 0;
702 crossscore[count+1][count+1] = 10000000.0;
703 cut1[count+1] = len1;
704 cut2[count+1] = len2;
705 count += 2;
706 count0 = count;
707
708 blockAlign2( cut1, cut2, sortedseg1, sortedseg2, crossscore, &count );
709 if( count0 > count )
710 {
711 #if 0
712 fprintf( stderr, "\7 REPEAT!? \n" );
713 #else
714 fprintf( stderr, "REPEAT!? \n" );
715 #endif
716 if( fftRepeatStop ) exit( 1 );
717 }
718 #if KEIKA
719 else fprintf( stderr, "done\n" );
720 #endif
721 }
722
723 #if 0
724 fftfp = fopen( "fft", "a" );
725 fprintf( fftfp, "RESULT after sort:\n" );
726 for( l=0; l<count; l++ )
727 {
728 fprintf( fftfp, "cut[%d]=%d, ", l, segment1[l].center );
729 fprintf( fftfp, "%d\n", segment2[l].center );
730 }
731 fclose( fftfp );
732 #endif
733
734 #if 0
735 fftfp = fopen( "fft", "a" );
736 fprintf( fftfp, "RESULT after sort:\n" );
737 for( l=0; l<count; l++ )
738 {
739 fprintf( fftfp, "cut : %d %d\n", cut1[l], cut2[l] );
740 }
741 fclose( fftfp );
742 #endif
743
744 #if KEIKA
745 fprintf( trap_g, "Devided to %d segments\n", count-1 );
746 fprintf( trap_g, "%d %d forg\n", MIN( clus1, clus2 ), count-1 );
747 #endif
748
749 totallen = 0;
750 for( j=0; j<clus1; j++ ) result1[j][0] = 0;
751 for( j=0; j<clus2; j++ ) result2[j][0] = 0;
752 totalscore = 0.0;
753 *totalimpmatch = 0.0;
754 for( i=0; i<count-1; i++ )
755 {
756 #if DEBUG
757 fprintf( stderr, "DP %03d / %03d %4d to ", i+1, count-1, totallen );
758 #else
759 #if KEIKA
760 fprintf( stderr, "DP %03d / %03d\r", i+1, count-1 );
761 #endif
762 #endif
763
764 if( cut1[i] )
765 {
766 getkyokaigap( sgap1, seq1, cut1[i]-1, clus1 );
767 getkyokaigap( sgap2, seq2, cut2[i]-1, clus2 );
768 }
769 else
770 {
771 for( j=0; j<clus1; j++ ) sgap1[j] = 'o';
772 for( j=0; j<clus2; j++ ) sgap2[j] = 'o';
773 }
774 if( cut1[i+1] != len1 )
775 {
776 getkyokaigap( egap1, seq1, cut1[i+1], clus1 );
777 getkyokaigap( egap2, seq2, cut2[i+1], clus2 );
778 }
779 else
780 {
781 for( j=0; j<clus1; j++ ) egap1[j] = 'o';
782 for( j=0; j<clus2; j++ ) egap2[j] = 'o';
783 }
784
785 for( j=0; j<clus1; j++ )
786 {
787 strncpy( tmpres1[j], seq1[j]+cut1[i], cut1[i+1]-cut1[i] );
788 tmpres1[j][cut1[i+1]-cut1[i]] = 0;
789 }
790 if( kobetsubunkatsu ) commongappick_record( clus1, tmpres1, gapmap1 );
791 for( j=0; j<clus2; j++ )
792 {
793 strncpy( tmpres2[j], seq2[j]+cut2[i], cut2[i+1]-cut2[i] );
794 tmpres2[j][cut2[i+1]-cut2[i]] = 0;
795 }
796 if( kobetsubunkatsu ) commongappick_record( clus2, tmpres2, gapmap2 );
797
798 #if 0
799 fprintf( stderr, "count = %d\n", count );
800 fprintf( stderr, "### reg1 = %d-%d\n", cut1[i], cut1[i+1]-1 );
801 fprintf( stderr, "### reg2 = %d-%d\n", cut2[i], cut2[i+1]-1 );
802 #endif
803
804 switch( alg )
805 {
806 case( 'a' ):
807 totalscore += Aalign( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen );
808 break;
809 case( 'A' ):
810 if( scoringmatrices ) // called by tditeration.c
811 {
812 totalscore += partA__align_variousdist( whichmtx, scoringmatrices, NULL, tmpres1, tmpres2, eff1, eff2, eff1s, eff2s, clus1, clus2, alloclen, localhom, &impmatch, cut1[i], cut1[i+1]-1, cut2[i], cut2[i+1]-1, gapmap1, gapmap2, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres );
813 }
814 else
815 totalscore += partA__align( tmpres1, tmpres2, eff1, eff2, clus1, clus2, alloclen, localhom, &impmatch, cut1[i], cut1[i+1]-1, cut2[i], cut2[i+1]-1, gapmap1, gapmap2, sgap1, sgap2, egap1, egap2, chudanpt, chudanref, chudanres );
816 *totalimpmatch += impmatch;
817 // fprintf( stderr, "*totalimpmatch in Falign_localhom = %f\n", *totalimpmatch );
818
819
820 break;
821 default:
822 fprintf( stderr, "alg = %c\n", alg );
823 ErrorExit( "ERROR IN SOURCE FILE Falign.c" );
824 break;
825 }
826 #ifdef enablemultithread
827 if( chudanres && *chudanres )
828 {
829 // fprintf( stderr, "\n\n## CHUUDAN!!! at Falign_localhom\n" );
830 return( -1.0 );
831 }
832 #endif
833
834 nlen = strlen( tmpres1[0] );
835 if( totallen + nlen > alloclen )
836 {
837 fprintf( stderr, "totallen=%d + nlen=%d > alloclen = %d\n", totallen, nlen, alloclen );
838 ErrorExit( "LENGTH OVER in Falign\n " );
839 }
840 for( j=0; j<clus1; j++ ) strcat( result1[j], tmpres1[j] );
841 for( j=0; j<clus2; j++ ) strcat( result2[j], tmpres2[j] );
842 totallen += nlen;
843 #if 0
844 fprintf( stderr, "%4d\r", totallen );
845 fprintf( stderr, "\n\n" );
846 for( j=0; j<clus1; j++ )
847 {
848 fprintf( stderr, "%s\n", tmpres1[j] );
849 }
850 fprintf( stderr, "-------\n" );
851 for( j=0; j<clus2; j++ )
852 {
853 fprintf( stderr, "%s\n", tmpres2[j] );
854 }
855 #endif
856 }
857 #if KEIKA
858 fprintf( stderr, "DP ... done \n" );
859 #endif
860
861 for( j=0; j<clus1; j++ ) strcpy( seq1[j], result1[j] );
862 for( j=0; j<clus2; j++ ) strcpy( seq2[j], result2[j] );
863 #if 0
864 for( j=0; j<clus1; j++ )
865 {
866 fprintf( stderr, "%s\n", result1[j] );
867 }
868 fprintf( stderr, "- - - - - - - - - - -\n" );
869 for( j=0; j<clus2; j++ )
870 {
871 fprintf( stderr, "%s\n", result2[j] );
872 }
873 #endif
874 return( totalscore );
875 }
876