Context Navigation

Falign_localhom.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: 21.1 KB

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

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

Context Navigation

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

Download in other formats: