3 #define SEGMENTSIZE 150
8 void keika( char *str, int current, int all )
11 fprintf( stderr, "%s : ", str );
13 fprintf( stderr, "\b\b\b\b\b\b\b\b" );
14 fprintf( stderr, "%3d /%3d", current+1, all+1 );
16 if( current+1 == all )
17 fprintf( stderr, "\b\b\b\b\b\b\b\bdone. \n" );
20 double maxItch( double *soukan, int size )
25 for( i=0; i<size; i++ )
26 if( ( cand = soukan[i] ) > value ) value = cand;
30 void calcNaiseki( Fukusosuu *value, Fukusosuu *x, Fukusosuu *y )
32 value->R = x->R * y->R + x->I * y->I;
33 value->I = -x->R * y->I + x->I * y->R;
36 Fukusosuu *AllocateFukusosuuVec( int l1 )
39 value = (Fukusosuu *)calloc( l1, sizeof( Fukusosuu ) );
42 fprintf( stderr, "Cannot allocate %d FukusosuuVec\n", l1 );
48 Fukusosuu **AllocateFukusosuuMtx( int l1, int l2 )
52 // fprintf( stderr, "allocating %d x %d FukusosuuMtx\n", l1, l2 );
53 value = (Fukusosuu **)calloc( l1+1, sizeof( Fukusosuu * ) );
56 fprintf( stderr, "Cannot allocate %d x %d FukusosuuVecMtx\n", l1, l2 );
61 value[j] = AllocateFukusosuuVec( l2 );
64 fprintf( stderr, "Cannot allocate %d x %d FukusosuuVecMtx\n", l1, l2 );
72 Fukusosuu ***AllocateFukusosuuCub( int l1, int l2, int l3 )
76 value = calloc( l1+1, sizeof( Fukusosuu ** ) );
77 if( !value ) ErrorExit( "Cannot allocate Fukusosuu" );
78 for( i=0; i<l1; i++ ) value[i] = AllocateFukusosuuMtx( l2, l3 );
83 void FreeFukusosuuVec( Fukusosuu *vec )
88 void FreeFukusosuuMtx( Fukusosuu **mtx )
92 for( i=0; mtx[i]; i++ )
93 free( (void *)mtx[i] );
97 int getKouho( int *kouho, int nkouho, double *soukan, int nlen2 )
100 int nlen4 = nlen2 / 2;
103 int ikouho = 0; // by D.Mathog, iinoka?
104 for( j=0; j<nkouho; j++ )
107 for( i=0; i<nlen2; i++ )
109 if( ( tmp = soukan[i] ) > max )
122 fprintf( stderr, "Kouho No.%d, pos=%d, score=%f, lag=%d\n", j, ikouho, soukan[ikouho], ikouho-nlen4 );
124 soukan[ikouho] = -9999.9;
125 kouho[j] = ( ikouho - nlen4 );
130 void zurasu2( int lag, int clus1, int clus2,
131 char **seq1, char **seq2,
132 char **aseq1, char **aseq2 )
136 fprintf( stderr, "### lag = %d\n", lag );
140 for( i=0; i<clus1; i++ ) aseq1[i] = seq1[i];
141 for( i=0; i<clus2; i++ ) aseq2[i] = seq2[i]+lag;
145 for( i=0; i<clus1; i++ ) aseq1[i] = seq1[i]-lag;
146 for( i=0; i<clus2; i++ ) aseq2[i] = seq2[i];
150 void zurasu( int lag, int clus1, int clus2,
151 char **seq1, char **seq2,
152 char **aseq1, char **aseq2 )
156 fprintf( stderr, "lag = %d\n", lag );
160 for( i=0; i<clus1; i++ ) strcpy( aseq1[i], seq1[i] );
161 for( i=0; i<clus2; i++ ) strcpy( aseq2[i], seq2[i]+lag );
165 for( i=0; i<clus1; i++ ) strcpy( aseq1[i], seq1[i]-lag );
166 for( i=0; i<clus2; i++ ) strcpy( aseq2[i], seq2[i] );
171 int alignableReagion( int clus1, int clus2,
172 char **seq1, char **seq2,
173 double *eff1, double *eff2,
177 int status, starttmp = 0; // by D.Mathog, a gess
181 int length = 0; // by D.Mathog, a gess
182 static TLS double *stra = NULL;
183 static TLS int alloclen = 0;
186 static TLS double threshold;
187 static TLS double *prf1 = NULL;
188 static TLS double *prf2 = NULL;
189 static TLS int *hat1 = NULL;
190 static TLS int *hat2 = NULL;
200 fprintf( stderr, "### In alignableRegion, clus1=%d, clus2=%d \n", clus1, clus2 );
201 fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
202 fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
203 fprintf( stderr, "eff1[0] = %f\n", eff1[0] );
204 fprintf( stderr, "eff2[0] = %f\n", eff2[0] );
209 FreeDoubleVec( stra ); stra = NULL;
210 FreeDoubleVec( prf1 ); prf1 = NULL;
211 FreeDoubleVec( prf2 ); prf2 = NULL;
212 FreeIntVec( hat1 ); hat1 = NULL;
213 FreeIntVec( hat2 ); hat2 = NULL;
219 prf1 = AllocateDoubleVec( 26 );
220 prf2 = AllocateDoubleVec( 26 );
221 hat1 = AllocateIntVec( 27 );
222 hat2 = AllocateIntVec( 27 );
225 len = MIN( strlen( seq1[0] ), strlen( seq2[0] ) );
226 maxlen = MAX( strlen( seq1[0] ), strlen( seq2[0] ) ) + fftWinSize;
227 if( alloclen < maxlen )
231 FreeDoubleVec( stra );
235 threshold = (int)fftThreshold / 100.0 * 600.0 * fftWinSize;
237 stra = AllocateDoubleVec( maxlen );
243 for( i=0; i<clus1; i++ ) for( j=0; j<clus2; j++ ) totaleff += eff1[i] * eff2[j];
244 for( i=0; i<len; i++ )
247 for( j=0; j<26; j++ )
256 while( j-- ) prf1[amino_n[(*seq1pt++)[i]]] += *eff1pt++;
258 for( j=0; j<clus1; j++ ) prf1[amino_n[(int)seq1[j][i]]] += eff1[j];
260 for( j=0; j<clus2; j++ ) prf2[amino_n[(int)seq2[j][i]]] += eff2[j];
264 for( j=25; j>=0; j-- )
280 /* make site score */
282 for( k=hat1[26]; k!=-1; k=hat1[k] )
283 for( j=hat2[26]; j!=-1; j=hat2[j] )
284 // stra[i] += n_dis[k][j] * prf1[k] * prf2[j];
285 stra[i] += n_disFFT[k][j] * prf1[k] * prf2[j];
289 (seg+0)->skipForeward = 0;
290 (seg+1)->skipBackward = 0;
294 for( j=0; j<fftWinSize; j++ ) score += stra[j];
296 for( i=1; i<len-fftWinSize; i++ )
298 score = score - stra[i-1] + stra[i+fftWinSize-1];
300 fprintf( stderr, "%d %10.0f ? %10.0f\n", i, score, threshold );
303 if( score > threshold )
308 seg->center = ( seg->start + seg->end + fftWinSize ) / 2 ;
324 if( score <= threshold || length > SEGMENTSIZE )
328 if( length > fftWinSize )
330 seg->start = starttmp;
332 seg->center = ( seg->start + seg->end + fftWinSize ) / 2 ;
333 seg->score = cumscore;
335 fprintf( stderr, "%d-%d length = %d, score = %f, value = %d\n", seg->start, seg->end, length, cumscore, value );
337 if( length > SEGMENTSIZE )
339 (seg+0)->skipForeward = 1;
340 (seg+1)->skipBackward = 1;
344 (seg+0)->skipForeward = 0;
345 (seg+1)->skipBackward = 0;
354 if( value > MAXSEG - 3 ) ErrorExit( "TOO MANY SEGMENTS!");
358 if( status && length > fftWinSize )
361 seg->start = starttmp;
362 seg->center = ( starttmp + i + fftWinSize ) / 2 ;
363 seg->score = cumscore;
365 fprintf( stderr, "%d-%d length = %d\n", seg->start, seg->end, length );
372 // fprintf( stderr, "returning %d\n", value );
377 static int permit( Segment *seg1, Segment *seg2 )
380 if( seg1->end >= seg2->start ) return( 0 );
381 if( seg1->pair->end >= seg2->pair->start ) return( 0 );
385 void blockAlign2( int *cut1, int *cut2, Segment **seg1, Segment **seg2, double **ocrossscore, int *ncut )
387 int i, j, k, shift, cur1, cur2, count, klim;
388 static TLS int crossscoresize = 0;
389 static TLS int *result1 = NULL;
390 static TLS int *result2 = NULL;
391 static TLS int *ocut1 = NULL;
392 static TLS int *ocut2 = NULL;
394 static TLS double **crossscore = NULL;
395 static TLS int **track = NULL;
396 static TLS double maxj, maxi;
397 static TLS int pointj, pointi;
408 FreeDoubleMtx( crossscore );
413 if( result1 == NULL )
415 result1 = AllocateIntVec( MAXSEG );
416 result2 = AllocateIntVec( MAXSEG );
417 ocut1 = AllocateIntVec( MAXSEG );
418 ocut2 = AllocateIntVec( MAXSEG );
421 if( crossscoresize < *ncut+2 )
423 crossscoresize = *ncut+2;
424 if( fftkeika ) fprintf( stderr, "allocating crossscore and track, size = %d\n", crossscoresize );
425 if( track ) FreeIntMtx( track );
426 if( crossscore ) FreeDoubleMtx( crossscore );
427 track = AllocateIntMtx( crossscoresize, crossscoresize );
428 crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
432 for( i=0; i<*ncut-2; i++ )
433 fprintf( stderr, "%d.start = %d, score = %f\n", i, seg1[i]->start, seg1[i]->score );
435 for( i=0; i<*ncut; i++ )
436 fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
437 for( i=0; i<*ncut; i++ )
439 for( j=0; j<*ncut; j++ )
440 fprintf( stderr, "%#4.0f ", ocrossscore[i][j] );
441 fprintf( stderr, "\n" );
445 for( i=0; i<*ncut; i++ ) for( j=0; j<*ncut; j++ ) /* mudadanaa */
446 crossscore[i][j] = ocrossscore[i][j];
447 for( i=0; i<*ncut; i++ )
453 for( i=1; i<*ncut; i++ )
456 fprintf( stderr, "### i=%d/%d\n", i,*ncut );
458 for( j=1; j<*ncut; j++ )
460 pointi = 0; maxi = 0.0;
462 for( k=0; k<klim; k++ )
465 fprintf( stderr, "k=%d, i=%d\n", k, i );
467 if( k && k<*ncut-1 && j<*ncut-1 && !permit( seg1[k-1], seg1[j-1] ) ) continue;
468 if( crossscore[i-1][k] > maxj )
471 maxi = crossscore[i-1][k];
475 pointj = 0; maxj = 0.0;
477 for( k=0; k<klim; k++ )
479 if( k && k<*ncut-1 && i<*ncut-1 && !permit( seg2[k-1], seg2[i-1] ) ) continue;
480 if( crossscore[k][j-1] > maxj )
483 maxj = crossscore[k][j-1];
490 maximum = crossscore[i-1][j-1];
496 track[i][j] = j - pointi;
502 track[i][j] = pointj - i;
505 crossscore[i][j] += maximum;
509 for( i=0; i<*ncut; i++ )
511 for( j=0; j<*ncut; j++ )
512 fprintf( stderr, "%3d ", track[i][j] );
513 fprintf( stderr, "\n" );
518 result1[MAXSEG-1] = *ncut-1;
519 result2[MAXSEG-1] = *ncut-1;
521 for( i=MAXSEG-1; i>=1; i-- )
525 if( cur1 == 0 || cur2 == 0 ) break;
526 shift = track[cur1][cur2];
529 result1[i-1] = cur1 - 1;
530 result2[i-1] = cur2 - 1;
535 result1[i-1] = cur1 - 1;
536 result2[i-1] = cur2 - shift;
540 result1[i-1] = cur1 + shift;
541 result2[i-1] = cur2 - 1;
546 for( j=i; j<MAXSEG; j++ )
548 if( ocrossscore[result1[j]][result2[j]] == 0.0 ) continue;
550 if( result1[j] == result1[j-1] || result2[j] == result2[j-1] )
551 if( ocrossscore[result1[j]][result2[j]] > ocrossscore[result1[j-1]][result2[j-1]] )
554 cut1[count] = ocut1[result1[j]];
555 cut2[count] = ocut2[result2[j]];
562 for( i=0; i<*ncut; i++ )
563 fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
567 void blockAlign3( int *cut1, int *cut2, Segment **seg1, Segment **seg2, double **ocrossscore, int *ncut )
568 // memory complexity = O(n^3), time complexity = O(n^2)
570 int i, j, shift, cur1, cur2, count;
571 static TLS int crossscoresize = 0;
572 static TLS int jumpposi, *jumppos;
573 static TLS double jumpscorei, *jumpscore;
574 static TLS int *result1 = NULL;
575 static TLS int *result2 = NULL;
576 static TLS int *ocut1 = NULL;
577 static TLS int *ocut2 = NULL;
579 static TLS double **crossscore = NULL;
580 static TLS int **track = NULL;
582 if( result1 == NULL )
584 result1 = AllocateIntVec( MAXSEG );
585 result2 = AllocateIntVec( MAXSEG );
586 ocut1 = AllocateIntVec( MAXSEG );
587 ocut2 = AllocateIntVec( MAXSEG );
589 if( crossscoresize < *ncut+2 )
591 crossscoresize = *ncut+2;
592 if( fftkeika ) fprintf( stderr, "allocating crossscore and track, size = %d\n", crossscoresize );
593 if( track ) FreeIntMtx( track );
594 if( crossscore ) FreeDoubleMtx( crossscore );
595 if( jumppos ) FreeIntVec( jumppos );
596 if( jumpscore ) FreeDoubleVec( jumpscore );
597 track = AllocateIntMtx( crossscoresize, crossscoresize );
598 crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
599 jumppos = AllocateIntVec( crossscoresize );
600 jumpscore = AllocateDoubleVec( crossscoresize );
604 for( i=0; i<*ncut-2; i++ )
605 fprintf( stderr, "%d.start = %d, score = %f\n", i, seg1[i]->start, seg1[i]->score );
607 for( i=0; i<*ncut; i++ )
608 fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
609 for( i=0; i<*ncut; i++ )
611 for( j=0; j<*ncut; j++ )
612 fprintf( stderr, "%#4.0f ", ocrossscore[i][j] );
613 fprintf( stderr, "\n" );
617 for( i=0; i<*ncut; i++ ) for( j=0; j<*ncut; j++ ) /* mudadanaa */
618 crossscore[i][j] = ocrossscore[i][j];
619 for( i=0; i<*ncut; i++ )
624 for( j=0; j<*ncut; j++ )
626 jumpscore[j] = -999.999;
630 for( i=1; i<*ncut; i++ )
633 jumpscorei = -999.999;
636 for( j=1; j<*ncut; j++ )
639 fprintf( stderr, "in blockalign3, ### i=%d, j=%d\n", i, j );
644 for( k=0; k<j-2; k++ )
647 fprintf( stderr, "k=%d, i=%d\n", k, i );
649 if( k && k<*ncut-1 && j<*ncut-1 && !permit( seg1[k-1], seg1[j-1] ) ) continue;
650 if( crossscore[i-1][k] > maxj )
653 maxi = crossscore[i-1][k];
657 pointj = 0; maxj = 0.0;
658 for( k=0; k<i-2; k++ )
660 if( k && k<*ncut-1 && i<*ncut-1 && !permit( seg2[k-1], seg2[i-1] ) ) continue;
661 if( crossscore[k][j-1] > maxj )
664 maxj = crossscore[k][j-1];
672 maximum = crossscore[i-1][j-1];
675 if( maximum < jumpscorei && permit( seg1[jumpposi], seg1[i] ) )
677 maximum = jumpscorei;
678 track[i][j] = j - jumpposi;
681 if( maximum < jumpscore[j] && permit( seg2[jumppos[j]], seg2[j] ) )
683 maximum = jumpscore[j];
684 track[i][j] = jumpscore[j] - i;
687 crossscore[i][j] += maximum;
689 if( jumpscorei < crossscore[i-1][j] )
691 jumpscorei = crossscore[i-1][j];
695 if( jumpscore[j] < crossscore[i][j-1] )
697 jumpscore[j] = crossscore[i][j-1];
703 for( i=0; i<*ncut; i++ )
705 for( j=0; j<*ncut; j++ )
706 fprintf( stderr, "%3d ", track[i][j] );
707 fprintf( stderr, "\n" );
712 result1[MAXSEG-1] = *ncut-1;
713 result2[MAXSEG-1] = *ncut-1;
715 for( i=MAXSEG-1; i>=1; i-- )
719 if( cur1 == 0 || cur2 == 0 ) break;
720 shift = track[cur1][cur2];
723 result1[i-1] = cur1 - 1;
724 result2[i-1] = cur2 - 1;
729 result1[i-1] = cur1 - 1;
730 result2[i-1] = cur2 - shift;
734 result1[i-1] = cur1 + shift;
735 result2[i-1] = cur2 - 1;
740 for( j=i; j<MAXSEG; j++ )
742 if( ocrossscore[result1[j]][result2[j]] == 0.0 ) continue;
744 if( result1[j] == result1[j-1] || result2[j] == result2[j-1] )
745 if( ocrossscore[result1[j]][result2[j]] > ocrossscore[result1[j-1]][result2[j-1]] )
748 cut1[count] = ocut1[result1[j]];
749 cut2[count] = ocut2[result2[j]];
756 for( i=0; i<*ncut; i++ )
757 fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );