Context Navigation

Falign.c

Visit:

Last change on this file was 10371, checked in by aboeckma, 12 years ago
updated mafft version. Added extensions (no svn ignore, yet)
File size: 60.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: branches/lib/GDE/MAFFT/mafft-7.055-with-extensions/core/Falign.c

Download in other formats: