Context Navigation

pairlara.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.3 KB

Line
1	#include "mltaln.h"
2
3	#define DEBUG 0
4	#define IODEBUG 0
5	#define SCOREOUT 0
6
7	static char *whereispairalign;
8	static char *laraparams;
9	static char foldalignopt[1000];
10
11	static void t2u( char *seq )
12	{
13	while( *seq )
14	{
15	if ( seq == 'A' ) seq = 'a';
16	else if( seq == 'a' ) seq = 'a';
17	else if( seq == 'T' ) seq = 'u';
18	else if( seq == 't' ) seq = 'u';
19	else if( seq == 'U' ) seq = 'u';
20	else if( seq == 'u' ) seq = 'u';
21	else if( seq == 'G' ) seq = 'g';
22	else if( seq == 'g' ) seq = 'g';
23	else if( seq == 'C' ) seq = 'c';
24	else if( seq == 'c' ) seq = 'c';
25	else *seq = 'n';
26	seq++;
27	}
28	}
29
30	static float recallpairfoldalign( char mseq1, char mseq2, int m1, int m2, int of1pt, int of2pt, int alloclen )
31	{
32	static FILE *fp = NULL;
33	float value;
34	char *aln1;
35	char *aln2;
36	int of1tmp, of2tmp;
37
38	if( fp == NULL )
39	{
40	fp = fopen( "_foldalignout", "r" );
41	if( fp == NULL )
42	{
43	fprintf( stderr, "Cannot open _foldalignout\n" );
44	exit( 1 );
45	}
46	}
47
48	aln1 = calloc( alloclen, sizeof( char ) );
49	aln2 = calloc( alloclen, sizeof( char ) );
50
51	readpairfoldalign( fp, mseq1, mseq2, aln1, aln2, m1, m2, &of1tmp, &of2tmp, alloclen );
52
53	if( strstr( foldalignopt, "-global") )
54	{
55	fprintf( stderr, "Calling G__align11\n" );
56	value = G__align11( mseq1, mseq2, alloclen );
57	*of1pt = 0;
58	*of2pt = 0;
59	}
60	else
61	{
62	fprintf( stderr, "Calling L__align11\n" );
63	value = L__align11( mseq1, mseq2, alloclen, of1pt, of2pt );
64	}
65
66	// value = (float)naivepairscore11( mseq1, mseq2, penalty ); // nennnotame
67
68	if( aln1[0] == 0 )
69	{
70	fprintf( stderr, "FOLDALIGN returned no alignment between %d and %d. Sequence alignment is used instead.\n", m1+1, m2+1 );
71	}
72	else
73	{
74	strcpy( *mseq1, aln1 );
75	strcpy( *mseq2, aln2 );
76	*of1pt = of1tmp;
77	*of2pt = of2tmp;
78	}
79
80	// value = naivepairscore11( mseq1, mseq2, penalty ); // v6.511 ha kore wo tsukau, global nomi dakara.
81
82	// fclose( fp ); // saigo dake yatta houga yoi.
83
84	// fprintf( stderr, "mseq1 = %s\n", mseq1 );
85	// fprintf( stderr, "mseq2 = %s\n", mseq2 );
86
87
88	free( aln1 );
89	free( aln2 );
90
91	return( value );
92	}
93
94	static void callfoldalign( int nseq, char **mseq )
95	{
96	FILE *fp;
97	int i;
98	int res;
99	static char com[10000];
100
101	for( i=0; i<nseq; i++ )
102	t2u( mseq[i] );
103
104	fp = fopen( "_foldalignin", "w" );
105	if( !fp )
106	{
107	fprintf( stderr, "Cannot open _foldalignin\n" );
108	exit( 1 );
109	}
110	for( i=0; i<nseq; i++ )
111	{
112	fprintf( fp, ">%d\n", i+1 );
113	fprintf( fp, "%s\n", mseq[i] );
114	}
115	fclose( fp );
116
117	sprintf( com, "env PATH=%s foldalign210 %s _foldalignin > _foldalignout ", whereispairalign, foldalignopt );
118	res = system( com );
119	if( res )
120	{
121	fprintf( stderr, "Error in foldalign\n" );
122	exit( 1 );
123	}
124
125	}
126
127	static void calllara( int nseq, char *mseq, char laraarg )
128	{
129	FILE *fp;
130	int i;
131	int res;
132	static char com[10000];
133
134	for( i=0; i<nseq; i++ )
135
136	fp = fopen( "_larain", "w" );
137	if( !fp )
138	{
139	fprintf( stderr, "Cannot open _larain\n" );
140	exit( 1 );
141	}
142	for( i=0; i<nseq; i++ )
143	{
144	fprintf( fp, ">%d\n", i+1 );
145	fprintf( fp, "%s\n", mseq[i] );
146	}
147	fclose( fp );
148
149
150	// fprintf( stderr, "calling LaRA\n" );
151	sprintf( com, "env PATH=%s:/bin:/usr/bin mafft_lara -i _larain -w _laraout -o _lara.params %s", whereispairalign, laraarg );
152	res = system( com );
153	if( res )
154	{
155	fprintf( stderr, "Error in lara\n" );
156	exit( 1 );
157	}
158	}
159
160	static float recalllara( char mseq1, char mseq2, int alloclen )
161	{
162	static FILE *fp = NULL;
163	static char *ungap1;
164	static char *ungap2;
165	static char *ori1;
166	static char *ori2;
167	int res;
168	static char com[10000];
169	float value;
170
171
172	if( fp == NULL )
173	{
174	fp = fopen( "_laraout", "r" );
175	if( fp == NULL )
176	{
177	fprintf( stderr, "Cannot open _laraout\n" );
178	exit( 1 );
179	}
180	ungap1 = AllocateCharVec( alloclen );
181	ungap2 = AllocateCharVec( alloclen );
182	ori1 = AllocateCharVec( alloclen );
183	ori2 = AllocateCharVec( alloclen );
184	}
185
186
187	strcpy( ori1, *mseq1 );
188	strcpy( ori2, *mseq2 );
189
190	fgets( com, 999, fp );
191	myfgets( com, 9999, fp );
192	strcpy( *mseq1, com );
193	myfgets( com, 9999, fp );
194	strcpy( *mseq2, com );
195
196	gappick0( ungap1, *mseq1 );
197	gappick0( ungap2, *mseq2 );
198	t2u( ungap1 );
199	t2u( ungap2 );
200
201	if( strcmp( ungap1, ori1 ) \|\| strcmp( ungap2, ori2 ) )
202	{
203	fprintf( stderr, "SEQUENCE CHANGED!!\n" );
204	fprintf( stderr, "mseq1 = %s\n", mseq1 );
205	fprintf( stderr, "ungap1 = %s\n", ungap1 );
206	fprintf( stderr, "ori1 = %s\n", ori1 );
207	fprintf( stderr, "mseq2 = %s\n", mseq2 );
208	fprintf( stderr, "ungap2 = %s\n", ungap2 );
209	fprintf( stderr, "ori2 = %s\n", ori2 );
210	exit( 1 );
211	}
212
213	value = (float)naivepairscore11( mseq1, mseq2, penalty );
214
215	// fclose( fp ); // saigo dake yatta houga yoi.
216
217	return( value );
218	}
219
220	static float callmxscarna( char mseq1, char mseq2, int alloclen )
221	{
222	FILE *fp;
223	int res;
224	static char com[10000];
225	float value;
226
227
228	t2u( *mseq1 );
229	t2u( *mseq2 );
230	fp = fopen( "_mxscarnain", "w" );
231	if( !fp )
232	{
233	fprintf( stderr, "Cannot open _mxscarnain\n" );
234	exit( 1 );
235	}
236	fprintf( fp, ">1\n" );
237	fprintf( fp, "%s\n", *mseq1 );
238	fprintf( fp, ">2\n" );
239	fprintf( fp, "%s\n", *mseq2 );
240	fclose( fp );
241
242	sprintf( com, "env PATH=%s mxscarna _mxscarnain > _mxscarnaout 2>/dev/null", whereispairalign );
243	res = system( com );
244	if( res )
245	{
246	fprintf( stderr, "Error in mxscarna\n" );
247	exit( 1 );
248	}
249
250	fp = fopen( "_mxscarnaout", "r" );
251	if( !fp )
252	{
253	fprintf( stderr, "Cannot open _mxscarnaout\n" );
254	exit( 1 );
255	}
256
257	fgets( com, 999, fp );
258	load1SeqWithoutName_new( fp, *mseq1 );
259	fgets( com, 999, fp );
260	load1SeqWithoutName_new( fp, *mseq2 );
261
262	fclose( fp );
263
264	// fprintf( stderr, "mseq1 = %s\n", mseq1 );
265	// fprintf( stderr, "mseq2 = %s\n", mseq2 );
266
267	value = (float)naivepairscore11( mseq1, mseq2, penalty );
268
269	return( value );
270	}
271
272
273	void arguments( int argc, char *argv[] )
274	{
275	int c;
276
277	foldalignopt[0] = 0;
278	laraparams = NULL;
279	inputfile = NULL;
280	fftkeika = 0;
281	pslocal = -1000.0;
282	constraint = 0;
283	nblosum = 62;
284	fmodel = 0;
285	calledByXced = 0;
286	devide = 0;
287	use_fft = 0;
288	fftscore = 1;
289	fftRepeatStop = 0;
290	fftNoAnchStop = 0;
291	weight = 3;
292	utree = 1;
293	tbutree = 1;
294	refine = 0;
295	check = 1;
296	cut = 0.0;
297	disp = 0;
298	outgap = 1;
299	alg = 'A';
300	mix = 0;
301	tbitr = 0;
302	scmtd = 5;
303	tbweight = 0;
304	tbrweight = 3;
305	checkC = 0;
306	treemethod = 'x';
307	contin = 0;
308	scoremtx = 1;
309	kobetsubunkatsu = 0;
310	divpairscore = 0;
311	dorp = NOTSPECIFIED;
312	ppenalty = NOTSPECIFIED;
313	ppenalty_OP = NOTSPECIFIED;
314	ppenalty_ex = NOTSPECIFIED;
315	ppenalty_EX = NOTSPECIFIED;
316	poffset = NOTSPECIFIED;
317	kimuraR = NOTSPECIFIED;
318	pamN = NOTSPECIFIED;
319	geta2 = GETA2;
320	fftWinSize = NOTSPECIFIED;
321	fftThreshold = NOTSPECIFIED;
322	RNAppenalty = NOTSPECIFIED;
323	RNApthr = NOTSPECIFIED;
324
325	while( --argc > 0 && (*++argv)[0] == '-' )
326	{
327	while ( ( c = *++argv[0] ) )
328	{
329	switch( c )
330	{
331	case 'i':
332	inputfile = *++argv;
333	fprintf( stderr, "inputfile = %s\n", inputfile );
334	--argc;
335	goto nextoption;
336	case 'f':
337	ppenalty = (int)( atof( ++argv ) 1000 - 0.5 );
338	--argc;
339	goto nextoption;
340	case 'g':
341	ppenalty_ex = (int)( atof( ++argv ) 1000 - 0.5 );
342	--argc;
343	goto nextoption;
344	case 'O':
345	ppenalty_OP = (int)( atof( ++argv ) 1000 - 0.5 );
346	--argc;
347	goto nextoption;
348	case 'E':
349	ppenalty_EX = (int)( atof( ++argv ) 1000 - 0.5 );
350	--argc;
351	goto nextoption;
352	case 'h':
353	poffset = (int)( atof( ++argv ) 1000 - 0.5 );
354	--argc;
355	goto nextoption;
356	case 'k':
357	kimuraR = atoi( *++argv );
358	// fprintf( stderr, "kimuraR = %d\n", kimuraR );
359	--argc;
360	goto nextoption;
361	case 'b':
362	nblosum = atoi( *++argv );
363	scoremtx = 1;
364	// fprintf( stderr, "blosum %d\n", nblosum );
365	--argc;
366	goto nextoption;
367	case 'j':
368	pamN = atoi( *++argv );
369	scoremtx = 0;
370	TMorJTT = JTT;
371	fprintf( stderr, "jtt %d\n", pamN );
372	--argc;
373	goto nextoption;
374	case 'm':
375	pamN = atoi( *++argv );
376	scoremtx = 0;
377	TMorJTT = TM;
378	fprintf( stderr, "TM %d\n", pamN );
379	--argc;
380	goto nextoption;
381	case 'l':
382	ppslocal = (int)( atof( ++argv ) 1000 + 0.5 );
383	pslocal = (int)( 600.0 / 1000.0 * ppslocal + 0.5);
384	// fprintf( stderr, "ppslocal = %d\n", ppslocal );
385	// fprintf( stderr, "pslocal = %d\n", pslocal );
386	--argc;
387	goto nextoption;
388	case 'd':
389	whereispairalign = *++argv;
390	fprintf( stderr, "whereispairalign = %s\n", whereispairalign );
391	--argc;
392	goto nextoption;
393	case 'p':
394	laraparams = *++argv;
395	fprintf( stderr, "laraparams = %s\n", laraparams );
396	--argc;
397	goto nextoption;
398	#if 1
399	case 'a':
400	fmodel = 1;
401	break;
402	#endif
403	case 'r':
404	fmodel = -1;
405	break;
406	case 'D':
407	dorp = 'd';
408	break;
409	case 'P':
410	dorp = 'p';
411	break;
412	case 'e':
413	fftscore = 0;
414	break;
415	#if 0
416	case 'O':
417	fftNoAnchStop = 1;
418	break;
419	#endif
420	case 'Q':
421	calledByXced = 1;
422	break;
423	case 'x':
424	disp = 1;
425	break;
426	#if 0
427	case 'a':
428	alg = 'a';
429	break;
430	#endif
431	case 'S':
432	alg = 'S';
433	break;
434	case 'L':
435	alg = 'L';
436	break;
437	case 's':
438	alg = 's';
439	break;
440	case 'B':
441	alg = 'B';
442	break;
443	case 'T':
444	alg = 'T';
445	break;
446	case 'H':
447	alg = 'H';
448	break;
449	case 'M':
450	alg = 'M';
451	break;
452	case 'R':
453	alg = 'R';
454	break;
455	case 'N':
456	alg = 'N';
457	break;
458	case 'A':
459	alg = 'A';
460	break;
461	case 'V':
462	alg = 'V';
463	break;
464	case 'C':
465	alg = 'C';
466	break;
467	case 'F':
468	use_fft = 1;
469	break;
470	case 'v':
471	tbrweight = 3;
472	break;
473	case 'y':
474	divpairscore = 1;
475	break;
476	/* Modified 01/08/27, default: user tree */
477	case 'J':
478	tbutree = 0;
479	break;
480	/* modification end. */
481	case 'o':
482	// foldalignopt = *++argv;
483	strcat( foldalignopt, " " );
484	strcat( foldalignopt, *++argv );
485	fprintf( stderr, "foldalignopt = %s\n", foldalignopt );
486	--argc;
487	goto nextoption;
488	case 'z':
489	fftThreshold = atoi( *++argv );
490	--argc;
491	goto nextoption;
492	case 'w':
493	fftWinSize = atoi( *++argv );
494	--argc;
495	goto nextoption;
496	case 'Z':
497	checkC = 1;
498	break;
499	default:
500	fprintf( stderr, "illegal option %c\n", c );
501	argc = 0;
502	break;
503	}
504	}
505	nextoption:
506	;
507	}
508	if( argc == 1 )
509	{
510	cut = atof( (*argv) );
511	argc--;
512	}
513	if( argc != 0 )
514	{
515	fprintf( stderr, "options: Check source file !\n" );
516	exit( 1 );
517	}
518	if( tbitr == 1 && outgap == 0 )
519	{
520	fprintf( stderr, "conflicting options : o, m or u\n" );
521	exit( 1 );
522	}
523	if( alg == 'C' && outgap == 0 )
524	{
525	fprintf( stderr, "conflicting options : C, o\n" );
526	exit( 1 );
527	}
528	}
529
530	int countamino( char *s, int end )
531	{
532	int val = 0;
533	while( end-- )
534	if( *s++ != '-' ) val++;
535	return( val );
536	}
537
538	static void pairalign( char name[M][B], int nlen[M], char seq, char aseq, char mseq1, char mseq2, double *effarr, int alloclen )
539	{
540	int i, j, ilim;
541	int clus1, clus2;
542	int off1, off2;
543	float pscore = 0.0; // by D.Mathog
544	static char indication1, indication2;
545	FILE hat2p, hat3p;
546	static double **distancemtx;
547	static double *effarr1 = NULL;
548	static double *effarr2 = NULL;
549	char *pt;
550	char *hat2file = "hat2";
551	LocalHom *localhomtable, tmpptr;
552	static char **pair;
553	int intdum;
554	double bunbo;
555
556	localhomtable = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
557	for( i=0; i<njob; i++)
558	{
559	localhomtable[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
560	for( j=0; j<njob; j++)
561	{
562	localhomtable[i][j].start1 = -1;
563	localhomtable[i][j].end1 = -1;
564	localhomtable[i][j].start2 = -1;
565	localhomtable[i][j].end2 = -1;
566	localhomtable[i][j].opt = -1.0;
567	localhomtable[i][j].next = NULL;
568	localhomtable[i][j].nokori = 0;
569	}
570	}
571
572	if( effarr1 == NULL )
573	{
574	distancemtx = AllocateDoubleMtx( njob, njob );
575	effarr1 = AllocateDoubleVec( njob );
576	effarr2 = AllocateDoubleVec( njob );
577	indication1 = AllocateCharVec( 150 );
578	indication2 = AllocateCharVec( 150 );
579	#if 0
580	#else
581	pair = AllocateCharMtx( njob, njob );
582	#endif
583	}
584
585	#if 0
586	fprintf( stderr, "##### fftwinsize = %d, fftthreshold = %d\n", fftWinSize, fftThreshold );
587	#endif
588
589	#if 0
590	for( i=0; i<njob; i++ )
591	fprintf( stderr, "TBFAST effarr[%d] = %f\n", i, effarr[i] );
592	#endif
593
594
595	// writePre( njob, name, nlen, aseq, 0 );
596
597	for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) pair[i][j] = 0;
598	for( i=0; i<njob; i++ ) pair[i][i] = 1;
599
600	if( alg == 'H' )
601	{
602	fprintf( stderr, "Calling FOLDALIGN with option '%s'\n", foldalignopt );
603	callfoldalign( njob, seq );
604	fprintf( stderr, "done.\n" );
605	}
606	if( alg == 'B' )
607	{
608	fprintf( stderr, "Calling LARA\n" );
609	calllara( njob, seq, "" );
610	fprintf( stderr, "done.\n" );
611	}
612	if( alg == 'T' )
613	{
614	fprintf( stderr, "Calling SLARA\n" );
615	calllara( njob, seq, "-s" );
616	fprintf( stderr, "done.\n" );
617	}
618
619	ilim = njob - 1;
620	for( i=0; i<ilim; i++ )
621	{
622	fprintf( stderr, "% 5d / %d\r", i, njob );
623	for( j=i+1; j<njob; j++ )
624	{
625
626	if( strlen( seq[i] ) == 0 \|\| strlen( seq[j] ) == 0 )
627	{
628	distancemtx[i][j] = pscore;
629	continue;
630	}
631
632	strcpy( aseq[i], seq[i] );
633	strcpy( aseq[j], seq[j] );
634	clus1 = conjuctionfortbfast_old( pair, i, aseq, mseq1, effarr1, effarr, indication1 );
635	clus2 = conjuctionfortbfast_old( pair, j, aseq, mseq2, effarr2, effarr, indication2 );
636	// fprintf( stderr, "mseq1 = %s\n", mseq1[0] );
637	// fprintf( stderr, "mseq2 = %s\n", mseq2[0] );
638
639	#if 0
640	fprintf( stderr, "group1 = %.66s", indication1 );
641	fprintf( stderr, "\n" );
642	fprintf( stderr, "group2 = %.66s", indication2 );
643	fprintf( stderr, "\n" );
644	#endif
645	// for( l=0; l<clus1; l++ ) fprintf( stderr, "## STEP-eff for mseq1-%d %f\n", l, effarr1[l] );
646
647	#if 1
648	if( use_fft )
649	{
650	pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum, NULL, 0, NULL );
651	off1 = off2 = 0;
652	}
653	else
654	#endif
655	{
656	switch( alg )
657	{
658	case( 'a' ):
659	pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
660	off1 = off2 = 0;
661	break;
662	case( 'A' ):
663	pscore = G__align11( mseq1, mseq2, alloclen, NULL, 0, NULL );
664	off1 = off2 = 0;
665	break;
666	#if 0
667	case( 'V' ):
668	pscore = VAalign11( mseq1, mseq2, alloclen, &off1, &off2, localhomtable[i]+j );
669	fprintf( stderr, "i,j = %d,%d, score = %f\n", i,j, pscore );
670	break;
671	case( 'S' ):
672	fprintf( stderr, "aligning %d-%d\n", i, j );
673	pscore = suboptalign11( mseq1, mseq2, alloclen, &off1, &off2, localhomtable[i]+j );
674	fprintf( stderr, "i,j = %d,%d, score = %f\n", i,j, pscore );
675	break;
676	#endif
677	case( 'N' ):
678	pscore = genL__align11( mseq1, mseq2, alloclen, &off1, &off2 );
679	// fprintf( stderr, "pscore = %f\n", pscore );
680	break;
681	case( 'L' ):
682	pscore = L__align11( mseq1, mseq2, alloclen, &off1, &off2 );
683	// fprintf( stderr, "pscore (1) = %f\n", pscore );
684	// pscore = (float)naivepairscore11( mseq1, mseq2, penalty ); // nennnotame
685	// fprintf( stderr, "pscore (2) = %f\n\n", pscore );
686	break;
687	case( 'H' ):
688	pscore = recallpairfoldalign( mseq1, mseq2, i, j, &off1, &off2, alloclen );
689	break;
690	case( 'B' ):
691	case( 'T' ):
692	pscore = recalllara( mseq1, mseq2, alloclen );
693	off1 = off2 = 0;
694	// fprintf( stderr, "lara, pscore = %f\n", pscore );
695	break;
696	case( 's' ):
697	pscore = callmxscarna( mseq1, mseq2, alloclen );
698	off1 = off2 = 0;
699	// fprintf( stderr, "scarna, pscore = %f\n", pscore );
700	break;
701	case( 'M' ):
702	// pscore = MSalign11( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL );
703	pscore = MSalign11( mseq1, mseq2, alloclen );
704	// fprintf( stderr, "pscore = %f\n", pscore );
705	break;
706	ErrorExit( "ERROR IN SOURCE FILE" );
707	}
708	}
709	distancemtx[i][j] = pscore;
710	#if SCOREOUT
711	fprintf( stderr, "score = %10.2f (%d,%d)\n", pscore, i, j );
712	#endif
713	// fprintf( stderr, "pslocal = %d\n", pslocal );
714	// offset = makelocal( mseq1, mseq2, pslocal );
715	#if 0
716	fprintf( stderr, "off1 = %d, off2 = %d\n", off1, off2 );
717	fprintf( stderr, ">%d\n%s\n>%d\n%s\n>\n", i, mseq1[0], j, mseq2[0] );
718	#endif
719
720	// putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j, countamino( mseq1, off1 ), countamino( mseq2, off2 ), pscore, strlen( mseq1[0] ) );
721	// fprintf( stderr, "pscore = %f\n", pscore );
722	if( alg == 'H' )
723	// if( alg == 'H' \|\| alg == 's' \|\| alg == 'B' ) // next version
724	putlocalhom_ext( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
725	else if( alg != 'S' && alg != 'V' )
726	putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
727	}
728	}
729	for( i=0; i<njob; i++ )
730	{
731	pscore = 0.0;
732	for( pt=seq[i]; *pt; pt++ )
733	pscore += amino_dis[(int)pt][(int)pt];
734	distancemtx[i][i] = pscore;
735
736	}
737
738	ilim = njob-1;
739	for( i=0; i<ilim; i++ )
740	{
741	for( j=i+1; j<njob; j++ )
742	{
743	bunbo = MIN( distancemtx[i][i], distancemtx[j][j] );
744	if( bunbo == 0.0 )
745	distancemtx[i][j] = 2.0;
746	else
747	distancemtx[i][j] = ( 1.0 - distancemtx[i][j] / bunbo ) * 2.0;
748	}
749	}
750
751	hat2p = fopen( hat2file, "w" );
752	if( !hat2p ) ErrorExit( "Cannot open hat2." );
753	WriteHat2( hat2p, njob, name, distancemtx );
754	fclose( hat2p );
755
756	fprintf( stderr, "##### writing hat3\n" );
757	hat3p = fopen( "hat3", "w" );
758	if( !hat3p ) ErrorExit( "Cannot open hat3." );
759	ilim = njob-1;
760	for( i=0; i<ilim; i++ )
761	{
762	for( j=i+1; j<njob; j++ )
763	{
764	for( tmpptr=localhomtable[i]+j; tmpptr; tmpptr=tmpptr->next )
765	{
766	if( tmpptr->opt == -1.0 ) continue;
767	fprintf( hat3p, "%d %d %d %7.5f %d %d %d %d %p\n", i, j, tmpptr->overlapaa, tmpptr->opt, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2, (void *)tmpptr->next );
768	}
769	}
770	}
771	fclose( hat3p );
772	#if DEBUG
773	fprintf( stderr, "calling FreeLocalHomTable\n" );
774	#endif
775	FreeLocalHomTable( localhomtable, njob );
776	#if DEBUG
777	fprintf( stderr, "done. FreeLocalHomTable\n" );
778	#endif
779	}
780
781	static void WriteOptions( FILE *fp )
782	{
783
784	if( dorp == 'd' ) fprintf( fp, "DNA\n" );
785	else
786	{
787	if ( scoremtx == 0 ) fprintf( fp, "JTT %dPAM\n", pamN );
788	else if( scoremtx == 1 ) fprintf( fp, "BLOSUM %d\n", nblosum );
789	else if( scoremtx == 2 ) fprintf( fp, "M-Y\n" );
790	}
791	fprintf( stderr, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
792	if( use_fft ) fprintf( fp, "FFT on\n" );
793
794	fprintf( fp, "tree-base method\n" );
795	if( tbrweight == 0 ) fprintf( fp, "unweighted\n" );
796	else if( tbrweight == 3 ) fprintf( fp, "clustalw-like weighting\n" );
797	if( tbitr \|\| tbweight )
798	{
799	fprintf( fp, "iterate at each step\n" );
800	if( tbitr && tbrweight == 0 ) fprintf( fp, " unweighted\n" );
801	if( tbitr && tbrweight == 3 ) fprintf( fp, " reversely weighted\n" );
802	if( tbweight ) fprintf( fp, " weighted\n" );
803	fprintf( fp, "\n" );
804	}
805
806	fprintf( fp, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
807
808	if( alg == 'a' )
809	fprintf( fp, "Algorithm A\n" );
810	else if( alg == 'A' )
811	fprintf( fp, "Algorithm A+\n" );
812	else if( alg == 'S' )
813	fprintf( fp, "Apgorithm S\n" );
814	else if( alg == 'C' )
815	fprintf( fp, "Apgorithm A+/C\n" );
816	else
817	fprintf( fp, "Unknown algorithm\n" );
818
819	if( use_fft )
820	{
821	fprintf( fp, "FFT on\n" );
822	if( dorp == 'd' )
823	fprintf( fp, "Basis : 4 nucleotides\n" );
824	else
825	{
826	if( fftscore )
827	fprintf( fp, "Basis : Polarity and Volume\n" );
828	else
829	fprintf( fp, "Basis : 20 amino acids\n" );
830	}
831	fprintf( fp, "Threshold of anchors = %d%%\n", fftThreshold );
832	fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
833	}
834	else
835	fprintf( fp, "FFT off\n" );
836	fflush( fp );
837	}
838
839
840	int main( int argc, char *argv[] )
841	{
842	static int nlen[M];
843	static char name[M][B], **seq;
844	static char mseq1, mseq2;
845	static char **aseq;
846	static char **bseq;
847	static double *eff;
848	int i;
849	FILE *infp;
850	char c;
851	int alloclen;
852
853	arguments( argc, argv );
854
855	if( inputfile )
856	{
857	infp = fopen( inputfile, "r" );
858	if( !infp )
859	{
860	fprintf( stderr, "Cannot open %s\n", inputfile );
861	exit( 1 );
862	}
863	}
864	else
865	infp = stdin;
866
867	getnumlen( infp );
868	rewind( infp );
869
870	if( njob < 2 )
871	{
872	fprintf( stderr, "At least 2 sequences should be input!\n"
873	"Only %d sequence found.\n", njob );
874	exit( 1 );
875	}
876	if( njob > M )
877	{
878	fprintf( stderr, "The number of sequences must be < %d\n", M );
879	fprintf( stderr, "Please try the splittbfast program for such large data.\n" );
880	exit( 1 );
881	}
882
883	seq = AllocateCharMtx( njob, nlenmax*9+1 );
884	aseq = AllocateCharMtx( njob, nlenmax*9+1 );
885	bseq = AllocateCharMtx( njob, nlenmax*9+1 );
886	mseq1 = AllocateCharMtx( njob, 0 );
887	mseq2 = AllocateCharMtx( njob, 0 );
888	alloclen = nlenmax*9;
889
890	eff = AllocateDoubleVec( njob );
891
892	#if 0
893	Read( name, nlen, seq );
894	#else
895	readData( infp, name, nlen, seq );
896	#endif
897	fclose( infp );
898
899	constants( njob, seq );
900
901	#if 0
902	fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
903	#endif
904
905	initSignalSM();
906
907	initFiles();
908
909	WriteOptions( trap_g );
910
911	c = seqcheck( seq );
912	if( c )
913	{
914	fprintf( stderr, "Illegal character %c\n", c );
915	exit( 1 );
916	}
917
918	// writePre( njob, name, nlen, seq, 0 );
919
920	for( i=0; i<njob; i++ ) eff[i] = 1.0;
921
922
923	for( i=0; i<njob; i++ ) gappick0( bseq[i], seq[i] );
924
925	pairalign( name, nlen, bseq, aseq, mseq1, mseq2, eff, alloclen );
926
927	fprintf( trap_g, "done.\n" );
928	#if DEBUG
929	fprintf( stderr, "closing trap_g\n" );
930	#endif
931	fclose( trap_g );
932
933	// writePre( njob, name, nlen, aseq, !contin );
934	#if 0
935	writeData( stdout, njob, name, nlen, aseq );
936	#endif
937	#if IODEBUG
938	fprintf( stderr, "OSHIMAI\n" );
939	#endif
940	SHOWVERSION;
941	return( 0 );
942	}

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

Context Navigation

source: trunk/GDE/MAFFT/mafft-7.055-with-extensions/core/pairlara.c

Download in other formats: