Context Navigation

pairlocalalign.c

Visit:

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

Line
1	#include "mltaln.h"
2
3	#define DEBUG 0
4	#define IODEBUG 0
5	#define SCOREOUT 0
6
7
8	#define NODIST -9999
9
10	static char *whereispairalign;
11	static char *laraparams;
12	static char foldalignopt[1000];
13	static int stdout_align;
14	static int stdout_dist;
15	static int store_localhom;
16	static int store_dist;
17	static int nadd;
18	static int laste;
19	static int lastm;
20	static int lastsubopt;
21	static int lastonce;
22
23	typedef struct _lastres
24	{
25	int score;
26	int start1;
27	int start2;
28	char *aln1;
29	char *aln2;
30	} Lastres;
31
32	typedef struct _reg
33	{
34	int start;
35	int end;
36	} Reg;
37
38	typedef struct _aln
39	{
40	int nreg;
41	Reg *reg1;
42	Reg *reg2;
43	} Aln;
44
45	typedef struct _lastresx
46	{
47	int score;
48	int naln;
49	Aln *aln;
50	} Lastresx;
51
52	#ifdef enablemultithread
53	typedef struct _jobtable
54	{
55	int i;
56	int j;
57	} Jobtable;
58
59	typedef struct _thread_arg
60	{
61	int thread_no;
62	int njob;
63	Jobtable *jobpospt;
64	char **name;
65	char **seq;
66	char **dseq;
67	int *thereisxineachseq;
68	LocalHom **localhomtable;
69	double **distancemtx;
70	double *selfscore;
71	char ***bpp;
72	Lastresx **lastresx;
73	int alloclen;
74	pthread_mutex_t *mutex_counter;
75	pthread_mutex_t *mutex_stdout;
76	} thread_arg_t;
77	#endif
78
79	typedef struct _lastcallthread_arg
80	{
81	int nq, nd;
82	char **dseq;
83	char **qseq;
84	Lastresx **lastresx;
85	#ifdef enablemultithread
86	int thread_no;
87	int *kshare;
88	pthread_mutex_t *mutex;
89	#endif
90	} lastcallthread_arg_t;
91
92	static void t2u( char *seq )
93	{
94	while( *seq )
95	{
96	if ( seq == 'A' ) seq = 'a';
97	else if( seq == 'a' ) seq = 'a';
98	else if( seq == 'T' ) seq = 'u';
99	else if( seq == 't' ) seq = 'u';
100	else if( seq == 'U' ) seq = 'u';
101	else if( seq == 'u' ) seq = 'u';
102	else if( seq == 'G' ) seq = 'g';
103	else if( seq == 'g' ) seq = 'g';
104	else if( seq == 'C' ) seq = 'c';
105	else if( seq == 'c' ) seq = 'c';
106	else *seq = 'n';
107	seq++;
108	}
109	}
110
111	static int removex( char d, char m )
112	{
113	int val = 0;
114	while( *m != 0 )
115	{
116	if( m == 'X' \|\| m == 'x' )
117	{
118	m++;
119	val++;
120	}
121	else
122	{
123	d++ = m++;
124	}
125	}
126	*d = 0;
127	return( val );
128	}
129
130	static void putlocalhom_last( char s1, char s2, LocalHom localhompt, Lastresx lastresx )
131	{
132	char pt1, pt2;
133	int naln, nreg;
134	int iscore;
135	int isumscore;
136	int sumoverlap;
137	LocalHom *tmppt = localhompt;
138	LocalHom *tmppt2;
139	LocalHom *localhompt0;
140	Reg rpt1, rpt2;
141	Aln *apt;
142	int nlocalhom = 0;
143	int len;
144
145	// fprintf( stderr, "s1=%s\n", s1 );
146	// fprintf( stderr, "s2=%s\n", s2 );
147
148
149	naln = lastresx->naln;
150	apt = lastresx->aln;
151
152	if( naln == 0 ) return;
153	while( naln-- )
154	{
155	rpt1 = apt->reg1;
156	rpt2 = apt->reg2;
157	nreg = apt->nreg;
158	isumscore = 0;
159	sumoverlap = 0;
160	while( nreg-- )
161	{
162	if( nlocalhom++ > 0 )
163	{
164	// fprintf( stderr, "reallocating ...\n" );
165	tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
166	// fprintf( stderr, "done\n" );
167	tmppt = tmppt->next;
168	tmppt->next = NULL;
169	}
170	tmppt->start1 = rpt1->start;
171	tmppt->start2 = rpt2->start;
172	tmppt->end1 = rpt1->end;
173	tmppt->end2 = rpt2->end;
174	if( rpt1 == apt->reg1 ) localhompt0 = tmppt; // ?
175
176	// fprintf( stderr, "in putlocalhom, reg1: %d-%d (nreg=%d)\n", rpt1->start, rpt1->end, lastresx->nreg );
177	// fprintf( stderr, "in putlocalhom, reg2: %d-%d (nreg=%d)\n", rpt2->start, rpt2->end, lastresx->nreg );
178
179	len = tmppt->end1 - tmppt->start1 + 1;
180
181	// fprintf( stderr, "tmppt->start1=%d\n", tmppt->start1 );
182	// fprintf( stderr, "tmppt->start2=%d\n", tmppt->start2 );
183
184	// fprintf( stderr, "s1+tmppt->start1=%.s\n", len, len, s1+tmppt->start1 );
185	// fprintf( stderr, "s2+tmppt->start2=%.s\n", len, len, s2+tmppt->start2 );
186
187	pt1 = s1 + tmppt->start1;
188	pt2 = s2 + tmppt->start2;
189	iscore = 0;
190	while( len-- )
191	{
192	iscore += n_dis[(int)amino_n[(int)pt1++]][(int)amino_n[(int)pt2++]]; // - offset $B$O$$$i$J$$$+$b(B
193	// fprintf( stderr, "len=%d, %c-%c, iscore(0) = %d\n", len, (pt1-1), (pt2-1), iscore );
194	}
195
196	if( divpairscore )
197	{
198	tmppt->overlapaa = tmppt->end2-tmppt->start2+1;
199	tmppt->opt = (double)iscore / tmppt->overlapaa * 5.8 / 600;
200	}
201	else
202	{
203	isumscore += iscore;
204	sumoverlap += tmppt->end2-tmppt->start2+1;
205	}
206	rpt1++;
207	rpt2++;
208	}
209	#if 0
210	fprintf( stderr, "iscore (1)= %d\n", iscore );
211	fprintf( stderr, "al1: %d - %d\n", start1, end1 );
212	fprintf( stderr, "al2: %d - %d\n", start2, end2 );
213	#endif
214
215	if( !divpairscore )
216	{
217	for( tmppt2=localhompt0; tmppt2; tmppt2=tmppt2->next )
218	{
219	tmppt2->overlapaa = sumoverlap;
220	tmppt2->opt = (double)isumscore * 5.8 / ( 600 * sumoverlap );
221	// fprintf( stderr, "tmpptr->opt = %f\n", tmppt->opt );
222	}
223	}
224	apt++;
225	}
226	}
227
228	static int countcomma( char *s )
229	{
230	int v = 0;
231	while( s ) if( s++ == ',' ) v++;
232	return( v );
233	}
234
235	static float recallpairfoldalign( char mseq1, char mseq2, int m1, int m2, int of1pt, int of2pt, int alloclen )
236	{
237	static FILE *fp = NULL;
238	float value;
239	char *aln1;
240	char *aln2;
241	int of1tmp, of2tmp;
242
243	if( fp == NULL )
244	{
245	fp = fopen( "_foldalignout", "r" );
246	if( fp == NULL )
247	{
248	fprintf( stderr, "Cannot open _foldalignout\n" );
249	exit( 1 );
250	}
251	}
252
253	aln1 = calloc( alloclen, sizeof( char ) );
254	aln2 = calloc( alloclen, sizeof( char ) );
255
256	readpairfoldalign( fp, mseq1, mseq2, aln1, aln2, m1, m2, &of1tmp, &of2tmp, alloclen );
257
258	if( strstr( foldalignopt, "-global") )
259	{
260	fprintf( stderr, "Calling G__align11\n" );
261	value = G__align11( mseq1, mseq2, alloclen, outgap, outgap );
262	*of1pt = 0;
263	*of2pt = 0;
264	}
265	else
266	{
267	fprintf( stderr, "Calling L__align11\n" );
268	value = L__align11( mseq1, mseq2, alloclen, of1pt, of2pt );
269	}
270
271	// value = (float)naivepairscore11( mseq1, mseq2, penalty ); // nennnotame
272
273	if( aln1[0] == 0 )
274	{
275	fprintf( stderr, "FOLDALIGN returned no alignment between %d and %d. Sequence alignment is used instead.\n", m1+1, m2+1 );
276	}
277	else
278	{
279	strcpy( *mseq1, aln1 );
280	strcpy( *mseq2, aln2 );
281	*of1pt = of1tmp;
282	*of2pt = of2tmp;
283	}
284
285	// value = naivepairscore11( mseq1, mseq2, penalty ); // v6.511 ha kore wo tsukau, global nomi dakara.
286
287	// fclose( fp ); // saigo dake yatta houga yoi.
288
289	// fprintf( stderr, "mseq1 = %s\n", mseq1 );
290	// fprintf( stderr, "mseq2 = %s\n", mseq2 );
291
292
293	free( aln1 );
294	free( aln2 );
295
296	return( value );
297	}
298
299	static void block2reg( char block, Reg reg1, Reg *reg2, int start1, int start2 )
300	{
301	Reg rpt1, rpt2;
302	char tpt, npt;
303	int pos1, pos2;
304	int len, glen1, glen2;
305	pos1 = start1;
306	pos2 = start2;
307	rpt1 = reg1;
308	rpt2 = reg2;
309	while( block )
310	{
311	block++;
312	// fprintf( stderr, "block = %s\n", block );
313	tpt = strchr( block, ':' );
314	npt = strchr( block, ',' );
315	if( !tpt \|\| tpt > npt )
316	{
317	len = atoi( block );
318	reg1->start = pos1;
319	reg2->start = pos2;
320	pos1 += len - 1;
321	pos2 += len - 1;
322	reg1->end = pos1;
323	reg2->end = pos2;
324	// fprintf( stderr, "in loop reg1: %d-%d\n", reg1->start, reg1->end );
325	// fprintf( stderr, "in loop reg2: %d-%d\n", reg2->start, reg2->end );
326	reg1++;
327	reg2++;
328	}
329	else
330	{
331	sscanf( block, "%d:%d", &glen1, &glen2 );
332	pos1 += glen1 + 1;
333	pos2 += glen2 + 1;
334	}
335	block = npt;
336
337	}
338	reg1->start = reg1->end = reg2->start = reg2->end = -1;
339
340	while( rpt1->start != -1 )
341	{
342	// fprintf( stderr, "reg1: %d-%d\n", rpt1->start, rpt1->end );
343	// fprintf( stderr, "reg2: %d-%d\n", rpt2->start, rpt2->end );
344	rpt1++;
345	rpt2++;
346	}
347	// apt1 = apt2 = 0;
348	// fprintf( stderr, "aln1 = %s\n", aln1 );
349	// fprintf( stderr, "aln2 = %s\n", aln2 );
350	}
351
352
353	static void readlastresx_singleq( FILE fp, int n1, int nameq, Lastresx *lastresx )
354	{
355	char *gett;
356	Aln *tmpaln;
357	int prevnaln, naln, nreg;
358	#if 0
359	int i, pstart, pend, end1, end2;
360	#endif
361	int score, name1, start1, alnSize1, seqSize1;
362	int name2, start2, alnSize2, seqSize2;
363	char strand1, strand2;
364	int includeintoscore;
365	gett = calloc( 10000, sizeof( char ) );
366
367	// fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
368	// fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
369
370	while( 1 )
371	{
372	fgets( gett, 9999, fp );
373	if( feof( fp ) ) break;
374	if( gett[0] == '#' ) continue;
375	// fprintf( stdout, "gett = %s\n", gett );
376	if( gett[strlen(gett)-1] != '\n' )
377	{
378	fprintf( stderr, "Too long line?\n" );
379	exit( 1 );
380	}
381
382	sscanf( gett, "%d %d %d %d %c %d %d %d %d %c %d",
383	&score, &name1, &start1, &alnSize1, &strand1, &seqSize1,
384	&name2, &start2, &alnSize2, &strand2, &seqSize2 );
385
386	if( alg == 'R' && name2 <= name1 ) continue;
387	if( name2 != nameq )
388	{
389	fprintf( stderr, "BUG!!!\n" );
390	exit( 1 );
391	}
392
393	// if( lastresx[name1][name2].score ) continue; // dame!!!!
394
395
396	prevnaln = lastresx[name1][name2].naln;
397	#if 0
398	for( i=0; i<prevnaln; i++ )
399	{
400	nreg = lastresx[name1][name2].aln[i].nreg;
401
402	pstart = lastresx[name1][name2].aln[i].reg1[0].start + 0;
403	pend = lastresx[name1][name2].aln[i].reg1[nreg-1].end - 0;
404	end1 = start1 + alnSize1;
405	// fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
406	if( pstart <= start1 && start1 <= pend && pend - start1 > 1 ) break;
407	if( pstart <= end1 && end1 <= pend && end1 - pstart > 1 ) break;
408
409	pstart = lastresx[name1][name2].aln[i].reg2[0].start + 0;
410	pend = lastresx[name1][name2].aln[i].reg2[nreg-1].end - 0;
411	end2 = start2 + alnSize2;
412	// fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
413	if( pstart <= start2 && start2 <= pend && pend - start2 > 1 ) break;
414	if( pstart <= end2 && end2 <= pend && end2 - pstart > 1 ) break;
415	}
416	includeintoscore = ( i == prevnaln );
417	#else
418	if( prevnaln ) includeintoscore = 0;
419	else includeintoscore = 1;
420	#endif
421	if( !includeintoscore && !lastsubopt )
422	continue;
423
424	naln = prevnaln + 1;
425	lastresx[name1][name2].naln = naln;
426	// fprintf( stderr, "OK! add this alignment to hat3, %d-%d, naln = %d->%d\n", name1, name2, prevnaln, naln );
427
428	if( ( tmpaln = (Aln )realloc( lastresx[name1][name2].aln, (naln) sizeof( Aln ) ) ) == NULL ) // yoyu nashi
429	{
430	fprintf( stderr, "Cannot reallocate lastresx[][].aln\n" );
431	exit( 1 );
432	}
433	else
434	lastresx[name1][name2].aln = tmpaln;
435
436	nreg = countcomma( gett )/2 + 1;
437	lastresx[name1][name2].aln[prevnaln].nreg = nreg;
438	// lastresx[name1][name2].aln[naln].nreg = -1;
439	// lastresx[name1][name2].aln[naln].reg1 = NULL;
440	// lastresx[name1][name2].aln[naln].reg2 = NULL;
441	// fprintf( stderr, "name1=%d, name2=%d, nreg=%d, prevnaln=%d\n", name1, name2, nreg, prevnaln );
442
443	if( ( lastresx[name1][name2].aln[prevnaln].reg1 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
444	{
445	fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
446	exit( 1 );
447	}
448
449	if( ( lastresx[name1][name2].aln[prevnaln].reg2 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
450	{
451	fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
452	exit( 1 );
453	}
454
455	// lastresx[name1][name2].aln[prevnaln].reg1[0].start = -1; // iranai?
456	// lastresx[name1][name2].aln[prevnaln].reg2[0].start = -1; // iranai?
457	block2reg( strrchr( gett, '\t' ), lastresx[name1][name2].aln[prevnaln].reg1, lastresx[name1][name2].aln[prevnaln].reg2, start1, start2 );
458
459	if( includeintoscore )
460	{
461	if( lastresx[name1][name2].score ) score += penalty;
462	lastresx[name1][name2].score += score;
463	}
464
465	// fprintf( stderr, "score(%d,%d) = %d\n", name1, name2, lastresx[name1][name2].score );
466	}
467	free( gett );
468	}
469
470	#ifdef enablemultithread
471	#if 0
472	static void readlastresx_group( FILE fp, Lastresx *lastresx )
473	{
474	char *gett;
475	Aln *tmpaln;
476	int prevnaln, naln, nreg;
477	#if 0
478	int i, pstart, pend, end1, end2;
479	#endif
480	int score, name1, start1, alnSize1, seqSize1;
481	int name2, start2, alnSize2, seqSize2;
482	char strand1, strand2;
483	int includeintoscore;
484	gett = calloc( 10000, sizeof( char ) );
485
486	// fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
487	// fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
488
489	while( 1 )
490	{
491	fgets( gett, 9999, fp );
492	if( feof( fp ) ) break;
493	if( gett[0] == '#' ) continue;
494	// fprintf( stdout, "gett = %s\n", gett );
495	if( gett[strlen(gett)-1] != '\n' )
496	{
497	fprintf( stderr, "Too long line?\n" );
498	exit( 1 );
499	}
500
501	sscanf( gett, "%d %d %d %d %c %d %d %d %d %c %d",
502	&score, &name1, &start1, &alnSize1, &strand1, &seqSize1,
503	&name2, &start2, &alnSize2, &strand2, &seqSize2 );
504
505	if( alg == 'R' && name2 <= name1 ) continue;
506
507	// if( lastresx[name1][name2].score ) continue; // dame!!!!
508
509	prevnaln = lastresx[name1][name2].naln;
510	#if 0
511	for( i=0; i<prevnaln; i++ )
512	{
513	nreg = lastresx[name1][name2].aln[i].nreg;
514
515	pstart = lastresx[name1][name2].aln[i].reg1[0].start + 0;
516	pend = lastresx[name1][name2].aln[i].reg1[nreg-1].end - 0;
517	end1 = start1 + alnSize1;
518	// fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
519	if( pstart <= start1 && start1 <= pend && pend - start1 > 3 ) break;
520	if( pstart <= end1 && end1 <= pend && end1 - pstart > 3 ) break;
521
522	pstart = lastresx[name1][name2].aln[i].reg2[0].start + 0;
523	pend = lastresx[name1][name2].aln[i].reg2[nreg-1].end - 0;
524	end2 = start2 + alnSize2;
525	// fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
526	if( pstart <= start2 && start2 <= pend && pend - start2 > 3 ) break;
527	if( pstart <= end2 && end2 <= pend && end2 - pstart > 3 ) break;
528	}
529	includeintoscore = ( i == prevnaln );
530	#else
531	if( prevnaln ) includeintoscore = 0;
532	else includeintoscore = 1;
533	#endif
534	if( !includeintoscore && !lastsubopt )
535	continue;
536
537	naln = prevnaln + 1;
538	lastresx[name1][name2].naln = naln;
539	// fprintf( stderr, "OK! add this alignment to hat3, %d-%d, naln = %d->%d\n", name1, name2, prevnaln, naln );
540
541	if( ( tmpaln = (Aln )realloc( lastresx[name1][name2].aln, (naln) sizeof( Aln ) ) ) == NULL ) // yoyu nashi
542	{
543	fprintf( stderr, "Cannot reallocate lastresx[][].aln\n" );
544	exit( 1 );
545	}
546	else
547	lastresx[name1][name2].aln = tmpaln;
548
549
550
551	nreg = countcomma( gett )/2 + 1;
552	lastresx[name1][name2].aln[prevnaln].nreg = nreg;
553	// lastresx[name1][name2].aln[naln].nreg = -1;
554	// lastresx[name1][name2].aln[naln].reg1 = NULL;
555	// lastresx[name1][name2].aln[naln].reg2 = NULL;
556	// fprintf( stderr, "name1=%d, name2=%d, nreg=%d, prevnaln=%d\n", name1, name2, nreg, prevnaln );
557
558	if( ( lastresx[name1][name2].aln[prevnaln].reg1 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
559	{
560	fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
561	exit( 1 );
562	}
563
564	if( ( lastresx[name1][name2].aln[prevnaln].reg2 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
565	{
566	fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
567	exit( 1 );
568	}
569
570	// lastresx[name1][name2].aln[prevnaln].reg1[0].start = -1; // iranai?
571	// lastresx[name1][name2].aln[prevnaln].reg2[0].start = -1; // iranai?
572	block2reg( strrchr( gett, '\t' ), lastresx[name1][name2].aln[prevnaln].reg1, lastresx[name1][name2].aln[prevnaln].reg2, start1, start2 );
573
574	if( includeintoscore )
575	{
576	if( lastresx[name1][name2].score ) score += penalty;
577	lastresx[name1][name2].score += score;
578	}
579
580	// fprintf( stderr, "score(%d,%d) = %d\n", name1, name2, lastresx[name1][name2].score );
581	}
582	free( gett );
583	}
584	#endif
585	#endif
586
587	static void readlastresx( FILE fp, int n1, int n2, Lastresx lastresx, char seq1, char *seq2 )
588	{
589	char *gett;
590	Aln *tmpaln;
591	int prevnaln, naln, nreg;
592	#if 0
593	int i, pstart, pend, end1, end2;
594	#endif
595	int score, name1, start1, alnSize1, seqSize1;
596	int name2, start2, alnSize2, seqSize2;
597	char strand1, strand2;
598	int includeintoscore;
599	gett = calloc( 10000, sizeof( char ) );
600
601	// fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
602	// fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
603
604	while( 1 )
605	{
606	fgets( gett, 9999, fp );
607	if( feof( fp ) ) break;
608	if( gett[0] == '#' ) continue;
609	// fprintf( stdout, "gett = %s\n", gett );
610	if( gett[strlen(gett)-1] != '\n' )
611	{
612	fprintf( stderr, "Too long line?\n" );
613	exit( 1 );
614	}
615
616	sscanf( gett, "%d %d %d %d %c %d %d %d %d %c %d",
617	&score, &name1, &start1, &alnSize1, &strand1, &seqSize1,
618	&name2, &start2, &alnSize2, &strand2, &seqSize2 );
619
620	if( alg == 'R' && name2 <= name1 ) continue;
621
622	// if( lastresx[name1][name2].score ) continue; // dame!!!!
623
624	prevnaln = lastresx[name1][name2].naln;
625	#if 0
626	for( i=0; i<prevnaln; i++ )
627	{
628	nreg = lastresx[name1][name2].aln[i].nreg;
629
630	pstart = lastresx[name1][name2].aln[i].reg1[0].start + 0;
631	pend = lastresx[name1][name2].aln[i].reg1[nreg-1].end - 0;
632	end1 = start1 + alnSize1;
633	// fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
634	if( pstart <= start1 && start1 <= pend && pend - start1 > 3 ) break;
635	if( pstart <= end1 && end1 <= pend && end1 - pstart > 3 ) break;
636
637	pstart = lastresx[name1][name2].aln[i].reg2[0].start + 0;
638	pend = lastresx[name1][name2].aln[i].reg2[nreg-1].end - 0;
639	end2 = start2 + alnSize2;
640	// fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
641	if( pstart <= start2 && start2 <= pend && pend - start2 > 3 ) break;
642	if( pstart <= end2 && end2 <= pend && end2 - pstart > 3 ) break;
643	}
644	includeintoscore = ( i == prevnaln );
645	#else
646	if( prevnaln ) includeintoscore = 0;
647	else includeintoscore = 1;
648	#endif
649	if( !includeintoscore && !lastsubopt )
650	continue;
651
652	naln = prevnaln + 1;
653	lastresx[name1][name2].naln = naln;
654	// fprintf( stderr, "OK! add this alignment to hat3, %d-%d, naln = %d->%d\n", name1, name2, prevnaln, naln );
655
656	if( ( tmpaln = (Aln )realloc( lastresx[name1][name2].aln, (naln) sizeof( Aln ) ) ) == NULL ) // yoyu nashi
657	{
658	fprintf( stderr, "Cannot reallocate lastresx[][].aln\n" );
659	exit( 1 );
660	}
661	else
662	lastresx[name1][name2].aln = tmpaln;
663
664
665
666	nreg = countcomma( gett )/2 + 1;
667	lastresx[name1][name2].aln[prevnaln].nreg = nreg;
668	// lastresx[name1][name2].aln[naln].nreg = -1;
669	// lastresx[name1][name2].aln[naln].reg1 = NULL;
670	// lastresx[name1][name2].aln[naln].reg2 = NULL;
671	// fprintf( stderr, "name1=%d, name2=%d, nreg=%d, prevnaln=%d\n", name1, name2, nreg, prevnaln );
672
673	if( ( lastresx[name1][name2].aln[prevnaln].reg1 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
674	{
675	fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
676	exit( 1 );
677	}
678
679	if( ( lastresx[name1][name2].aln[prevnaln].reg2 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
680	{
681	fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
682	exit( 1 );
683	}
684
685	// lastresx[name1][name2].aln[prevnaln].reg1[0].start = -1; // iranai?
686	// lastresx[name1][name2].aln[prevnaln].reg2[0].start = -1; // iranai?
687	block2reg( strrchr( gett, '\t' ), lastresx[name1][name2].aln[prevnaln].reg1, lastresx[name1][name2].aln[prevnaln].reg2, start1, start2 );
688
689	if( includeintoscore )
690	{
691	if( lastresx[name1][name2].score ) score += penalty;
692	lastresx[name1][name2].score += score;
693	}
694
695	// fprintf( stderr, "score(%d,%d) = %d\n", name1, name2, lastresx[name1][name2].score );
696	}
697	free( gett );
698	}
699
700	#ifdef enablemultithread
701	#if 0
702	static void lastcallthread_group( void arg )
703	{
704	lastcallthread_arg_t targ = (lastcallthread_arg_t )arg;
705	int k, i;
706	int nq = targ->nq;
707	int nd = targ->nd;
708	#ifdef enablemultithread
709	int thread_no = targ->thread_no;
710	int *kshare = targ->kshare;
711	#endif
712	Lastresx **lastresx = targ->lastresx;
713	char **dseq = targ->dseq;
714	char **qseq = targ->qseq;
715	char command[5000];
716	FILE *lfp;
717	int msize;
718	int klim;
719	int qstart, qend, shou, amari;
720	char kd[1000];
721
722	if( nthread )
723	{
724	shou = nq / nthread;
725	amari = nq - shou * nthread;
726	fprintf( stderr, "shou: %d, amari: %d\n", shou, amari );
727
728	qstart = thread_no * shou;
729	if( thread_no - 1 < amari ) qstart += thread_no;
730	else qstart += amari;
731
732	qend = qstart + shou - 1;
733	if( thread_no < amari ) qend += 1;
734	fprintf( stderr, "%d: %d-%d\n", thread_no, qstart, qend );
735	}
736	k = -1;
737	while( 1 )
738	{
739	if( nthread )
740	{
741	if( qstart > qend ) break;
742	if( k == thread_no ) break;
743	fprintf( stderr, "\n%d-%d / %d (thread %d) \n", qstart, qend, nq, thread_no );
744	k = thread_no;
745	}
746	else
747	{
748	k++;
749	if( k == nq ) break;
750	fprintf( stderr, "\r%d / %d \r", k, nq );
751	}
752
753	if( alg == 'R' ) // if 'r' -> calllast_fast
754	{
755	fprintf( stderr, "Not supported\n" );
756	exit( 1 );
757	}
758	else // 'r'
759	{
760	kd[0] = 0;
761	}
762
763	sprintf( command, "_q%d", k );
764	lfp = fopen( command, "w" );
765	if( !lfp )
766	{
767	fprintf( stderr, "Cannot open %s", command );
768	exit( 1 );
769	}
770	for( i=qstart; i<=qend; i++ )
771	fprintf( lfp, ">%d\n%s\n", i, qseq[i] );
772	fclose( lfp );
773
774	// if( alg == 'R' ) msize = MAX(10,k+nq);
775	// else msize = MAX(10,nd+nq);
776	if( alg == 'R' ) msize = MAX(10,k*lastm);
777	else msize = MAX(10,nd*lastm);
778
779	// fprintf( stderr, "Calling lastal from lastcallthread, msize = %d, k=%d\n", msize, k );
780	// sprintf( command, "grep '>' _db%sd", kd );
781	// system( command );
782	sprintf( command, "%s/lastal -m %d -e %d -f 0 -s 1 -p _scoringmatrixforlast -a %d -b %d _db%sd _q%d > _lastres%d", whereispairalign, msize, laste, -penalty, -penalty_ex, kd, k, k );
783	if( system( command ) ) exit( 1 );
784
785	sprintf( command, "_lastres%d", k );
786	lfp = fopen( command, "r" );
787	if( !lfp )
788	{
789	fprintf( stderr, "Cannot read _lastres%d", k );
790	exit( 1 );
791	}
792	// readlastres( lfp, nd, nq, lastres, dseq, qseq );
793	// fprintf( stderr, "Reading lastres\n" );
794	readlastresx_group( lfp, lastresx );
795	fclose( lfp );
796	}
797	return( NULL );
798	}
799	#endif
800	#endif
801
802	static void lastcallthread( void arg )
803	{
804	lastcallthread_arg_t targ = (lastcallthread_arg_t )arg;
805	int k, i;
806	int nq = targ->nq;
807	int nd = targ->nd;
808	#ifdef enablemultithread
809	int thread_no = targ->thread_no;
810	int *kshare = targ->kshare;
811	#endif
812	Lastresx **lastresx = targ->lastresx;
813	char **dseq = targ->dseq;
814	char **qseq = targ->qseq;
815	char command[5000];
816	FILE *lfp;
817	int msize;
818	int klim;
819	char kd[1000];
820
821	k = -1;
822	while( 1 )
823	{
824
825	#ifdef enablemultithread
826	if( nthread )
827	{
828	pthread_mutex_lock( targ->mutex );
829	k = *kshare;
830	if( k == nq )
831	{
832	pthread_mutex_unlock( targ->mutex );
833	break;
834	}
835	fprintf( stderr, "\r%d / %d (thread %d) \r", k, nq, thread_no );
836	++(*kshare);
837	pthread_mutex_unlock( targ->mutex );
838	}
839	else
840	#endif
841	{
842	k++;
843	if( k == nq ) break;
844	fprintf( stderr, "\r%d / %d \r", k, nq );
845	}
846
847	if( alg == 'R' ) // if 'r' -> calllast_fast
848	{
849	klim = MIN( k, njob-nadd );
850	// klim = k; // dochira demo yoi
851	if( klim == k )
852	{
853	sprintf( command, "_db%dd", k );
854	lfp = fopen( command, "w" );
855	if( !lfp )
856	{
857	fprintf( stderr, "Cannot open _db." );
858	exit( 1 );
859	}
860	for( i=0; i<klim; i++ ) fprintf( lfp, ">%d\n%s\n", i, dseq[i] );
861	fclose( lfp );
862
863	// sprintf( command, "md5sum _db%dd > /dev/tty", k );
864	// system( command );
865
866	if( dorp == 'd' )
867	sprintf( command, "%s/lastdb _db%dd _db%dd", whereispairalign, k, k );
868	else
869	sprintf( command, "%s/lastdb -p _db%dd _db%dd", whereispairalign, k, k );
870	system( command );
871	sprintf( kd, "%d", k );
872	}
873	else // calllast_fast de tsukutta nowo riyou
874	{
875	kd[0] = 0;
876	// fprintf( stderr, "klim=%d, njob=%d, nadd=%d, skip!\n", klim, njob, nadd );
877	}
878	}
879	else // 'r'
880	{
881	kd[0] = 0;
882	}
883
884	sprintf( command, "_q%d", k );
885	lfp = fopen( command, "w" );
886	if( !lfp )
887	{
888	fprintf( stderr, "Cannot open %s", command );
889	exit( 1 );
890	}
891	fprintf( lfp, ">%d\n%s\n", k, qseq[k] );
892	fclose( lfp );
893
894	// if( alg == 'R' ) msize = MAX(10,k+nq);
895	// else msize = MAX(10,nd+nq);
896	if( alg == 'R' ) msize = MAX(10,k*lastm);
897	else msize = MAX(10,nd*lastm);
898
899	// fprintf( stderr, "Calling lastal from lastcallthread, msize = %d, k=%d\n", msize, k );
900	// sprintf( command, "grep '>' _db%sd", kd );
901	// system( command );
902	sprintf( command, "%s/lastal -m %d -e %d -f 0 -s 1 -p _scoringmatrixforlast -a %d -b %d _db%sd _q%d > _lastres%d", whereispairalign, msize, laste, -penalty, -penalty_ex, kd, k, k );
903	if( system( command ) ) exit( 1 );
904
905	sprintf( command, "_lastres%d", k );
906	lfp = fopen( command, "r" );
907	if( !lfp )
908	{
909	fprintf( stderr, "Cannot read _lastres%d", k );
910	exit( 1 );
911	}
912	// readlastres( lfp, nd, nq, lastres, dseq, qseq );
913	// fprintf( stderr, "Reading lastres\n" );
914	readlastresx_singleq( lfp, nd, k, lastresx );
915	fclose( lfp );
916	}
917	return( NULL );
918	}
919
920
921	static void calllast_fast( int nd, char dseq, int nq, char qseq, Lastresx **lastresx )
922	{
923	int i, j;
924	FILE *lfp;
925	char command[1000];
926
927	lfp = fopen( "_scoringmatrixforlast", "w" );
928	if( !lfp )
929	{
930	fprintf( stderr, "Cannot open _scoringmatrixforlast" );
931	exit( 1 );
932	}
933	if( dorp == 'd' )
934	{
935	fprintf( lfp, " " );
936	for( j=0; j<4; j++ ) fprintf( lfp, " %c ", amino[j] );
937	fprintf( lfp, "\n" );
938	for( i=0; i<4; i++ )
939	{
940	fprintf( lfp, "%c ", amino[i] );
941	for( j=0; j<4; j++ ) fprintf( lfp, " %d ", n_dis[i][j] );
942	fprintf( lfp, "\n" );
943	}
944	}
945	else
946	{
947	fprintf( lfp, " " );
948	for( j=0; j<20; j++ ) fprintf( lfp, " %c ", amino[j] );
949	fprintf( lfp, "\n" );
950	for( i=0; i<20; i++ )
951	{
952	fprintf( lfp, "%c ", amino[i] );
953	for( j=0; j<20; j++ ) fprintf( lfp, " %d ", n_dis[i][j] );
954	fprintf( lfp, "\n" );
955	}
956	}
957	fclose( lfp );
958
959	// if( alg == 'r' ) // if 'R' -> lastcallthread, kokonoha nadd>0 no toki nomi shiyou
960	{
961	sprintf( command, "_dbd" );
962	lfp = fopen( command, "w" );
963	if( !lfp )
964	{
965	fprintf( stderr, "Cannot open _dbd" );
966	exit( 1 );
967	}
968	if( alg == 'R' )
969	j = njob-nadd;
970	else
971	j = nd;
972	for( i=0; i<j; i++ ) fprintf( lfp, ">%d\n%s\n", i, dseq[i] );
973
974	fclose( lfp );
975	if( dorp == 'd' )
976	sprintf( command, "%s/lastdb _dbd _dbd", whereispairalign );
977	else
978	sprintf( command, "%s/lastdb -p _dbd _dbd", whereispairalign );
979	system( command );
980	}
981
982	#ifdef enablemultithread
983	if( nthread )
984	{
985	pthread_t *handle;
986	pthread_mutex_t mutex;
987	lastcallthread_arg_t *targ;
988	int *ksharept;
989	targ = (lastcallthread_arg_t *)calloc( nthread, sizeof( lastcallthread_arg_t ) );
990	handle = calloc( nthread, sizeof( pthread_t ) );
991	ksharept = calloc( 1, sizeof(int) );
992	*ksharept = 0;
993	pthread_mutex_init( &mutex, NULL );
994	for( i=0; i<nthread; i++ )
995	{
996	targ[i].thread_no = i;
997	targ[i].kshare = ksharept;
998	targ[i].nq = nq;
999	targ[i].nd = nd;
1000	targ[i].dseq = dseq;
1001	targ[i].qseq = qseq;
1002	targ[i].lastresx = lastresx;
1003	targ[i].mutex = &mutex;
1004	pthread_create( handle+i, NULL, lastcallthread, (void *)(targ+i) );
1005	}
1006
1007	for( i=0; i<nthread; i++ )
1008	{
1009	pthread_join( handle[i], NULL );
1010	}
1011	pthread_mutex_destroy( &mutex );
1012	free( handle );
1013	free( targ );
1014	free( ksharept );
1015	}
1016	else
1017	#endif
1018	{
1019	lastcallthread_arg_t *targ;
1020	targ = (lastcallthread_arg_t *)calloc( 1, sizeof( lastcallthread_arg_t ) );
1021	targ[0].nq = nq;
1022	targ[0].nd = nd;
1023	targ[0].dseq = dseq;
1024	targ[0].qseq = qseq;
1025	targ[0].lastresx = lastresx;
1026	lastcallthread( targ );
1027	free( targ );
1028	}
1029
1030	}
1031
1032	static void calllast_once( int nd, char dseq, int nq, char qseq, Lastresx **lastresx )
1033	{
1034	int i, j;
1035	char command[5000];
1036	FILE *lfp;
1037	int msize;
1038	int res;
1039
1040	fprintf( stderr, "nq=%d\n", nq );
1041
1042	lfp = fopen( "_db", "w" );
1043	if( !lfp )
1044	{
1045	fprintf( stderr, "Cannot open _db" );
1046	exit( 1 );
1047	}
1048	for( i=0; i<nd; i++ ) fprintf( lfp, ">%d\n%s\n", i, dseq[i] );
1049	fclose( lfp );
1050
1051	if( dorp == 'd' )
1052	{
1053	sprintf( command, "%s/lastdb _db _db", whereispairalign );
1054	system( command );
1055	lfp = fopen( "_scoringmatrixforlast", "w" );
1056	if( !lfp )
1057	{
1058	fprintf( stderr, "Cannot open _scoringmatrixforlast" );
1059	exit( 1 );
1060	}
1061	fprintf( lfp, " " );
1062	for( j=0; j<4; j++ ) fprintf( lfp, " %c ", amino[j] );
1063	fprintf( lfp, "\n" );
1064	for( i=0; i<4; i++ )
1065	{
1066	fprintf( lfp, "%c ", amino[i] );
1067	for( j=0; j<4; j++ ) fprintf( lfp, " %d ", n_dis[i][j] );
1068	fprintf( lfp, "\n" );
1069	}
1070	fclose( lfp );
1071	#if 0
1072	sprintf( command, "lastex -s 2 -a %d -b %d -p _scoringmatrixforlast -E 10000 _db.prj _db.prj > _lastex", -penalty, -penalty_ex );
1073	system( command );
1074	lfp = fopen( "_lastex", "r" );
1075	fgets( command, 4999, lfp );
1076	fgets( command, 4999, lfp );
1077	fgets( command, 4999, lfp );
1078	fgets( command, 4999, lfp );
1079	laste = atoi( command );
1080	fclose( lfp );
1081	fprintf( stderr, "laste = %d\n", laste );
1082	sleep( 10 );
1083	#else
1084	// laste = 5000;
1085	#endif
1086	}
1087	else
1088	{
1089	sprintf( command, "%s/lastdb -p _db _db", whereispairalign );
1090	system( command );
1091	lfp = fopen( "_scoringmatrixforlast", "w" );
1092	if( !lfp )
1093	{
1094	fprintf( stderr, "Cannot open _scoringmatrixforlast" );
1095	exit( 1 );
1096	}
1097	fprintf( lfp, " " );
1098	for( j=0; j<20; j++ ) fprintf( lfp, " %c ", amino[j] );
1099	fprintf( lfp, "\n" );
1100	for( i=0; i<20; i++ )
1101	{
1102	fprintf( lfp, "%c ", amino[i] );
1103	for( j=0; j<20; j++ ) fprintf( lfp, " %d ", n_dis[i][j] );
1104	fprintf( lfp, "\n" );
1105	}
1106	fclose( lfp );
1107	// fprintf( stderr, "Not written yet\n" );
1108	}
1109
1110	lfp = fopen( "_q", "w" );
1111	if( !lfp )
1112	{
1113	fprintf( stderr, "Cannot open _q" );
1114	exit( 1 );
1115	}
1116	for( i=0; i<nq; i++ )
1117	{
1118	fprintf( lfp, ">%d\n%s\n", i, qseq[i] );
1119	}
1120	fclose( lfp );
1121
1122	msize = MAX(10,nd*lastm);
1123
1124	// fprintf( stderr, "Calling lastal from calllast_once, msize=%d\n", msize );
1125	sprintf( command, "%s/lastal -v -m %d -e %d -f 0 -s 1 -p _scoringmatrixforlast -a %d -b %d _db _q > _lastres", whereispairalign, msize, laste, -penalty, -penalty_ex );
1126	// sprintf( command, "lastal -v -m %d -e %d -f 0 -s 1 -p _scoringmatrixforlast -a %d -b %d _db _q > _lastres", 1, laste, -penalty, -penalty_ex );
1127	// sprintf( command, "lastal -v -e 40 -f 0 -s 1 -p _scoringmatrixforlast -a %d -b %d _db _q > _lastres", -penalty, -penalty_ex );
1128	res = system( command );
1129	if( res )
1130	{
1131	fprintf( stderr, "LAST aborted\n" );
1132	exit( 1 );
1133	}
1134
1135	lfp = fopen( "_lastres", "r" );
1136	if( !lfp )
1137	{
1138	fprintf( stderr, "Cannot read _lastres" );
1139	exit( 1 );
1140	}
1141	// readlastres( lfp, nd, nq, lastres, dseq, qseq );
1142	fprintf( stderr, "Reading lastres\n" );
1143	readlastresx( lfp, nd, nq, lastresx, dseq, qseq );
1144	fclose( lfp );
1145	}
1146
1147	static void callfoldalign( int nseq, char **mseq )
1148	{
1149	FILE *fp;
1150	int i;
1151	int res;
1152	static char com[10000];
1153
1154	for( i=0; i<nseq; i++ )
1155	t2u( mseq[i] );
1156
1157	fp = fopen( "_foldalignin", "w" );
1158	if( !fp )
1159	{
1160	fprintf( stderr, "Cannot open _foldalignin\n" );
1161	exit( 1 );
1162	}
1163	for( i=0; i<nseq; i++ )
1164	{
1165	fprintf( fp, ">%d\n", i+1 );
1166	fprintf( fp, "%s\n", mseq[i] );
1167	}
1168	fclose( fp );
1169
1170	sprintf( com, "env PATH=%s foldalign210 %s _foldalignin > _foldalignout ", whereispairalign, foldalignopt );
1171	res = system( com );
1172	if( res )
1173	{
1174	fprintf( stderr, "Error in foldalign\n" );
1175	exit( 1 );
1176	}
1177
1178	}
1179
1180	static void calllara( int nseq, char *mseq, char laraarg )
1181	{
1182	FILE *fp;
1183	int i;
1184	int res;
1185	static char com[10000];
1186
1187	// for( i=0; i<nseq; i++ )
1188
1189	fp = fopen( "_larain", "w" );
1190	if( !fp )
1191	{
1192	fprintf( stderr, "Cannot open _larain\n" );
1193	exit( 1 );
1194	}
1195	for( i=0; i<nseq; i++ )
1196	{
1197	fprintf( fp, ">%d\n", i+1 );
1198	fprintf( fp, "%s\n", mseq[i] );
1199	}
1200	fclose( fp );
1201
1202
1203	// fprintf( stderr, "calling LaRA\n" );
1204	sprintf( com, "env PATH=%s:/bin:/usr/bin mafft_lara -i _larain -w _laraout -o _lara.params %s", whereispairalign, laraarg );
1205	res = system( com );
1206	if( res )
1207	{
1208	fprintf( stderr, "Error in lara\n" );
1209	exit( 1 );
1210	}
1211	}
1212
1213	static float recalllara( char mseq1, char mseq2, int alloclen )
1214	{
1215	static FILE *fp = NULL;
1216	static char *ungap1;
1217	static char *ungap2;
1218	static char *ori1;
1219	static char *ori2;
1220	// int res;
1221	static char com[10000];
1222	float value;
1223
1224
1225	if( fp == NULL )
1226	{
1227	fp = fopen( "_laraout", "r" );
1228	if( fp == NULL )
1229	{
1230	fprintf( stderr, "Cannot open _laraout\n" );
1231	exit( 1 );
1232	}
1233	ungap1 = AllocateCharVec( alloclen );
1234	ungap2 = AllocateCharVec( alloclen );
1235	ori1 = AllocateCharVec( alloclen );
1236	ori2 = AllocateCharVec( alloclen );
1237	}
1238
1239
1240	strcpy( ori1, *mseq1 );
1241	strcpy( ori2, *mseq2 );
1242
1243	fgets( com, 999, fp );
1244	myfgets( com, 9999, fp );
1245	strcpy( *mseq1, com );
1246	myfgets( com, 9999, fp );
1247	strcpy( *mseq2, com );
1248
1249	gappick0( ungap1, *mseq1 );
1250	gappick0( ungap2, *mseq2 );
1251	t2u( ungap1 );
1252	t2u( ungap2 );
1253	t2u( ori1 );
1254	t2u( ori2 );
1255
1256	if( strcmp( ungap1, ori1 ) \|\| strcmp( ungap2, ori2 ) )
1257	{
1258	fprintf( stderr, "SEQUENCE CHANGED!!\n" );
1259	fprintf( stderr, "mseq1 = %s\n", mseq1 );
1260	fprintf( stderr, "ungap1 = %s\n", ungap1 );
1261	fprintf( stderr, "ori1 = %s\n", ori1 );
1262	fprintf( stderr, "mseq2 = %s\n", mseq2 );
1263	fprintf( stderr, "ungap2 = %s\n", ungap2 );
1264	fprintf( stderr, "ori2 = %s\n", ori2 );
1265	exit( 1 );
1266	}
1267
1268	value = (float)naivepairscore11( mseq1, mseq2, penalty );
1269
1270	// fclose( fp ); // saigo dake yatta houga yoi.
1271
1272	return( value );
1273	}
1274
1275
1276	static float calldafs_giving_bpp( char mseq1, char mseq2, char bpp1, char bpp2, int alloclen, int i, int j )
1277	{
1278	FILE *fp;
1279	int res;
1280	char *com;
1281	float value;
1282	char *dirname;
1283
1284
1285	dirname = calloc( 100, sizeof( char ) );
1286	com = calloc( 1000, sizeof( char ) );
1287	sprintf( dirname, "_%d-%d", i, j );
1288	sprintf( com, "rm -rf %s", dirname );
1289	system( com );
1290	sprintf( com, "mkdir %s", dirname );
1291	system( com );
1292
1293
1294	sprintf( com, "%s/_bpporg", dirname );
1295	fp = fopen( com, "w" );
1296	if( !fp )
1297	{
1298	fprintf( stderr, "Cannot write to %s/_bpporg\n", dirname );
1299	exit( 1 );
1300	}
1301	fprintf( fp, ">a\n" );
1302	while( *bpp1 )
1303	fprintf( fp, "%s", *bpp1++ );
1304
1305	fprintf( fp, ">b\n" );
1306	while( *bpp2 )
1307	fprintf( fp, "%s", *bpp2++ );
1308	fclose( fp );
1309
1310	sprintf( com, "tr -d '\\r' < %s/_bpporg > %s/_bpp", dirname, dirname );
1311	system( com ); // for cygwin, wakaran
1312
1313	t2u( *mseq1 );
1314	t2u( *mseq2 );
1315
1316	sprintf( com, "%s/_dafsinorg", dirname );
1317	fp = fopen( com, "w" );
1318	if( !fp )
1319	{
1320	fprintf( stderr, "Cannot open %s/_dafsinorg\n", dirname );
1321	exit( 1 );
1322	}
1323	fprintf( fp, ">1\n" );
1324	// fprintf( fp, "%s\n", *mseq1 );
1325	write1seq( fp, *mseq1 );
1326	fprintf( fp, ">2\n" );
1327	// fprintf( fp, "%s\n", *mseq2 );
1328	write1seq( fp, *mseq2 );
1329	fclose( fp );
1330
1331	sprintf( com, "tr -d '\\r' < %s/_dafsinorg > %s/_dafsin", dirname, dirname );
1332	system( com ); // for cygwin, wakaran
1333
1334	sprintf( com, "_dafssh%s", dirname );
1335	fp = fopen( com, "w" );
1336	fprintf( fp, "cd %s\n", dirname );
1337	fprintf( fp, "%s/dafs --mafft-in _bpp _dafsin > _dafsout 2>_dum\n", whereispairalign );
1338	fprintf( fp, "exit $tatus\n" );
1339	fclose( fp );
1340
1341	sprintf( com, "tr -d '\\r' < _dafssh%s > _dafssh%s.unix", dirname, dirname );
1342	system( com ); // for cygwin, wakaran
1343
1344	sprintf( com, "sh _dafssh%s.unix 2>_dum%s", dirname, dirname );
1345	res = system( com );
1346	if( res )
1347	{
1348	fprintf( stderr, "Error in dafs\n" );
1349	exit( 1 );
1350	}
1351
1352	sprintf( com, "%s/_dafsout", dirname );
1353
1354	fp = fopen( com, "r" );
1355	if( !fp )
1356	{
1357	fprintf( stderr, "Cannot open %s/_dafsout\n", dirname );
1358	exit( 1 );
1359	}
1360
1361	myfgets( com, 999, fp ); // nagai kanousei ga arunode
1362	fgets( com, 999, fp );
1363	myfgets( com, 999, fp ); // nagai kanousei ga arunode
1364	fgets( com, 999, fp );
1365	load1SeqWithoutName_new( fp, *mseq1 );
1366	fgets( com, 999, fp );
1367	load1SeqWithoutName_new( fp, *mseq2 );
1368
1369	fclose( fp );
1370
1371	// fprintf( stderr, "mseq1 = %s\n", mseq1 );
1372	// fprintf( stderr, "mseq2 = %s\n", mseq2 );
1373
1374	value = (float)naivepairscore11( mseq1, mseq2, penalty );
1375
1376	#if 0
1377	sprintf( com, "rm -rf %s > /dev/null 2>&1", dirname );
1378	if( system( com ) )
1379	{
1380	fprintf( stderr, "retrying to rmdir\n" );
1381	usleep( 2000 );
1382	system( com );
1383	}
1384	#endif
1385
1386	free( dirname );
1387	free( com );
1388
1389
1390	return( value );
1391	}
1392
1393	static float callmxscarna_giving_bpp( char mseq1, char mseq2, char bpp1, char bpp2, int alloclen, int i, int j )
1394	{
1395	FILE *fp;
1396	int res;
1397	char *com;
1398	float value;
1399	char *dirname;
1400
1401
1402	dirname = calloc( 100, sizeof( char ) );
1403	com = calloc( 1000, sizeof( char ) );
1404	sprintf( dirname, "_%d-%d", i, j );
1405	sprintf( com, "rm -rf %s", dirname );
1406	system( com );
1407	sprintf( com, "mkdir %s", dirname );
1408	system( com );
1409
1410
1411	sprintf( com, "%s/_bpporg", dirname );
1412	fp = fopen( com, "w" );
1413	if( !fp )
1414	{
1415	fprintf( stderr, "Cannot write to %s/_bpporg\n", dirname );
1416	exit( 1 );
1417	}
1418	fprintf( fp, ">a\n" );
1419	while( *bpp1 )
1420	fprintf( fp, "%s", *bpp1++ );
1421
1422	fprintf( fp, ">b\n" );
1423	while( *bpp2 )
1424	fprintf( fp, "%s", *bpp2++ );
1425	fclose( fp );
1426
1427	sprintf( com, "tr -d '\\r' < %s/_bpporg > %s/_bpp", dirname, dirname );
1428	system( com ); // for cygwin, wakaran
1429
1430	t2u( *mseq1 );
1431	t2u( *mseq2 );
1432
1433	sprintf( com, "%s/_mxscarnainorg", dirname );
1434	fp = fopen( com, "w" );
1435	if( !fp )
1436	{
1437	fprintf( stderr, "Cannot open %s/_mxscarnainorg\n", dirname );
1438	exit( 1 );
1439	}
1440	fprintf( fp, ">1\n" );
1441	// fprintf( fp, "%s\n", *mseq1 );
1442	write1seq( fp, *mseq1 );
1443	fprintf( fp, ">2\n" );
1444	// fprintf( fp, "%s\n", *mseq2 );
1445	write1seq( fp, *mseq2 );
1446	fclose( fp );
1447
1448	sprintf( com, "tr -d '\\r' < %s/_mxscarnainorg > %s/_mxscarnain", dirname, dirname );
1449	system( com ); // for cygwin, wakaran
1450
1451	#if 0
1452	sprintf( com, "cd %s; %s/mxscarnamod -readbpp _mxscarnain > _mxscarnaout 2>_dum", dirname, whereispairalign );
1453	#else
1454	sprintf( com, "_mxscarnash%s", dirname );
1455	fp = fopen( com, "w" );
1456	fprintf( fp, "cd %s\n", dirname );
1457	fprintf( fp, "%s/mxscarnamod -readbpp _mxscarnain > _mxscarnaout 2>_dum\n", whereispairalign );
1458	fprintf( fp, "exit $tatus\n" );
1459	fclose( fp );
1460	//sleep( 10000 );
1461
1462	sprintf( com, "tr -d '\\r' < _mxscarnash%s > _mxscarnash%s.unix", dirname, dirname );
1463	system( com ); // for cygwin, wakaran
1464
1465	sprintf( com, "sh _mxscarnash%s.unix 2>_dum%s", dirname, dirname );
1466	#endif
1467	res = system( com );
1468	if( res )
1469	{
1470	fprintf( stderr, "Error in mxscarna\n" );
1471	exit( 1 );
1472	}
1473
1474	sprintf( com, "%s/_mxscarnaout", dirname );
1475
1476	fp = fopen( com, "r" );
1477	if( !fp )
1478	{
1479	fprintf( stderr, "Cannot open %s/_mxscarnaout\n", dirname );
1480	exit( 1 );
1481	}
1482
1483	fgets( com, 999, fp );
1484	load1SeqWithoutName_new( fp, *mseq1 );
1485	fgets( com, 999, fp );
1486	load1SeqWithoutName_new( fp, *mseq2 );
1487
1488	fclose( fp );
1489
1490	// fprintf( stderr, "mseq1 = %s\n", mseq1 );
1491	// fprintf( stderr, "mseq2 = %s\n", mseq2 );
1492
1493	value = (float)naivepairscore11( mseq1, mseq2, penalty );
1494
1495	#if 0
1496	sprintf( com, "rm -rf %s > /dev/null 2>&1", dirname );
1497	if( system( com ) )
1498	{
1499	fprintf( stderr, "retrying to rmdir\n" );
1500	usleep( 2000 );
1501	system( com );
1502	}
1503	#endif
1504
1505	free( dirname );
1506	free( com );
1507
1508
1509	return( value );
1510	}
1511
1512	static void readhat4( FILE fp, char **bpp )
1513	{
1514	char oneline[1000];
1515	int bppsize;
1516	int onechar;
1517	// double prob;
1518	// int posi, posj;
1519
1520	bppsize = 0;
1521	// fprintf( stderr, "reading hat4\n" );
1522	onechar = getc(fp);
1523	// fprintf( stderr, "onechar = %c\n", onechar );
1524	if( onechar != '>' )
1525	{
1526	fprintf( stderr, "Format error\n" );
1527	exit( 1 );
1528	}
1529	ungetc( onechar, fp );
1530	fgets( oneline, 999, fp );
1531	while( 1 )
1532	{
1533	onechar = getc(fp);
1534	ungetc( onechar, fp );
1535	if( onechar == '>' \|\| onechar == EOF )
1536	{
1537	// fprintf( stderr, "Next\n" );
1538	bpp = realloc( bpp, (bppsize+2) * sizeof( char * ) );
1539	(*bpp)[bppsize] = NULL;
1540	break;
1541	}
1542	fgets( oneline, 999, fp );
1543	// fprintf( stderr, "oneline=%s\n", oneline );
1544	// sscanf( oneline, "%d %d %f", &posi, &posj, &prob );
1545	// fprintf( stderr, "%d %d -> %f\n", posi, posj, prob );
1546	bpp = realloc( bpp, (bppsize+2) * sizeof( char * ) );
1547	(*bpp)[bppsize] = calloc( 100, sizeof( char ) );
1548	strcpy( (*bpp)[bppsize], oneline );
1549	bppsize++;
1550	}
1551	}
1552
1553	static void preparebpp( int nseq, char ***bpp )
1554	{
1555	FILE *fp;
1556	int i;
1557
1558	fp = fopen( "hat4", "r" );
1559	if( !fp )
1560	{
1561	fprintf( stderr, "Cannot open hat4\n" );
1562	exit( 1 );
1563	}
1564	for( i=0; i<nseq; i++ )
1565	readhat4( fp, bpp+i );
1566	fclose( fp );
1567	}
1568
1569	void arguments( int argc, char *argv[] )
1570	{
1571	int c;
1572
1573	nthread = 1;
1574	laste = 5000;
1575	lastm = 3;
1576	nadd = 0;
1577	lastsubopt = 0;
1578	lastonce = 0;
1579	foldalignopt[0] = 0;
1580	laraparams = NULL;
1581	inputfile = NULL;
1582	fftkeika = 0;
1583	pslocal = -1000.0;
1584	constraint = 0;
1585	nblosum = 62;
1586	fmodel = 0;
1587	calledByXced = 0;
1588	devide = 0;
1589	use_fft = 0;
1590	fftscore = 1;
1591	fftRepeatStop = 0;
1592	fftNoAnchStop = 0;
1593	weight = 3;
1594	utree = 1;
1595	tbutree = 1;
1596	refine = 0;
1597	check = 1;
1598	cut = 0.0;
1599	disp = 0;
1600	outgap = 1;
1601	alg = 'A';
1602	mix = 0;
1603	tbitr = 0;
1604	scmtd = 5;
1605	tbweight = 0;
1606	tbrweight = 3;
1607	checkC = 0;
1608	treemethod = 'x';
1609	contin = 0;
1610	scoremtx = 1;
1611	kobetsubunkatsu = 0;
1612	divpairscore = 0;
1613	stdout_align = 0;
1614	stdout_dist = 0;
1615	store_dist = 1;
1616	store_localhom = 1;
1617	dorp = NOTSPECIFIED;
1618	ppenalty = NOTSPECIFIED;
1619	ppenalty_OP = NOTSPECIFIED;
1620	ppenalty_ex = NOTSPECIFIED;
1621	ppenalty_EX = NOTSPECIFIED;
1622	poffset = NOTSPECIFIED;
1623	kimuraR = NOTSPECIFIED;
1624	pamN = NOTSPECIFIED;
1625	geta2 = GETA2;
1626	fftWinSize = NOTSPECIFIED;
1627	fftThreshold = NOTSPECIFIED;
1628	RNAppenalty = NOTSPECIFIED;
1629	RNApthr = NOTSPECIFIED;
1630
1631	while( --argc > 0 && (*++argv)[0] == '-' )
1632	{
1633	while ( ( c = *++argv[0] ) )
1634	{
1635	switch( c )
1636	{
1637	case 'i':
1638	inputfile = *++argv;
1639	fprintf( stderr, "inputfile = %s\n", inputfile );
1640	--argc;
1641	goto nextoption;
1642	case 'f':
1643	ppenalty = (int)( atof( ++argv ) 1000 - 0.5 );
1644	--argc;
1645	goto nextoption;
1646	case 'g':
1647	ppenalty_ex = (int)( atof( ++argv ) 1000 - 0.5 );
1648	--argc;
1649	goto nextoption;
1650	case 'O':
1651	ppenalty_OP = (int)( atof( ++argv ) 1000 - 0.5 );
1652	--argc;
1653	goto nextoption;
1654	case 'E':
1655	ppenalty_EX = (int)( atof( ++argv ) 1000 - 0.5 );
1656	--argc;
1657	goto nextoption;
1658	case 'h':
1659	poffset = (int)( atof( ++argv ) 1000 - 0.5 );
1660	--argc;
1661	goto nextoption;
1662	case 'k':
1663	kimuraR = atoi( *++argv );
1664	// fprintf( stderr, "kimuraR = %d\n", kimuraR );
1665	--argc;
1666	goto nextoption;
1667	case 'b':
1668	nblosum = atoi( *++argv );
1669	scoremtx = 1;
1670	// fprintf( stderr, "blosum %d\n", nblosum );
1671	--argc;
1672	goto nextoption;
1673	case 'j':
1674	pamN = atoi( *++argv );
1675	scoremtx = 0;
1676	TMorJTT = JTT;
1677	fprintf( stderr, "jtt %d\n", pamN );
1678	--argc;
1679	goto nextoption;
1680	case 'm':
1681	pamN = atoi( *++argv );
1682	scoremtx = 0;
1683	TMorJTT = TM;
1684	fprintf( stderr, "TM %d\n", pamN );
1685	--argc;
1686	goto nextoption;
1687	#if 0
1688	case 'l':
1689	ppslocal = (int)( atof( ++argv ) 1000 + 0.5 );
1690	pslocal = (int)( 600.0 / 1000.0 * ppslocal + 0.5);
1691	// fprintf( stderr, "ppslocal = %d\n", ppslocal );
1692	// fprintf( stderr, "pslocal = %d\n", pslocal );
1693	--argc;
1694	goto nextoption;
1695	#else
1696	case 'l':
1697	if( atof( *++argv ) < 0.00001 ) store_localhom = 0;
1698	--argc;
1699	goto nextoption;
1700	#endif
1701	case 'd':
1702	whereispairalign = *++argv;
1703	fprintf( stderr, "whereispairalign = %s\n", whereispairalign );
1704	--argc;
1705	goto nextoption;
1706	case 'p':
1707	laraparams = *++argv;
1708	fprintf( stderr, "laraparams = %s\n", laraparams );
1709	--argc;
1710	goto nextoption;
1711	case 'C':
1712	nthread = atoi( *++argv );
1713	fprintf( stderr, "nthread = %d\n", nthread );
1714	--argc;
1715	goto nextoption;
1716	case 'I':
1717	nadd = atoi( *++argv );
1718	fprintf( stderr, "nadd = %d\n", nadd );
1719	--argc;
1720	goto nextoption;
1721	case 'w':
1722	lastm = atoi( *++argv );
1723	fprintf( stderr, "lastm = %d\n", lastm );
1724	--argc;
1725	goto nextoption;
1726	case 'e':
1727	laste = atoi( *++argv );
1728	fprintf( stderr, "laste = %d\n", laste );
1729	--argc;
1730	goto nextoption;
1731	case 'K': // Hontou ha iranai. disttbfast.c, tbfast.c to awaserutame.
1732	break;
1733	case 'c':
1734	stdout_dist = 1;
1735	break;
1736	case 'n':
1737	stdout_align = 1;
1738	break;
1739	case 'x':
1740	store_localhom = 0;
1741	store_dist = 0;
1742	break;
1743	#if 1
1744	case 'a':
1745	fmodel = 1;
1746	break;
1747	#endif
1748	#if 0
1749	case 'r':
1750	fmodel = -1;
1751	break;
1752	#endif
1753	case 'D':
1754	dorp = 'd';
1755	break;
1756	case 'P':
1757	dorp = 'p';
1758	break;
1759	#if 0
1760	case 'e':
1761	fftscore = 0;
1762	break;
1763	case 'O':
1764	fftNoAnchStop = 1;
1765	break;
1766	#endif
1767	#if 0
1768	case 'Q':
1769	calledByXced = 1;
1770	break;
1771	case 'x':
1772	disp = 1;
1773	break;
1774	case 'a':
1775	alg = 'a';
1776	break;
1777	case 'S':
1778	alg = 'S';
1779	break;
1780	#endif
1781	case 'Q':
1782	lastonce = 1;
1783	break;
1784	case 'S':
1785	lastsubopt = 1;
1786	break;
1787	case 't':
1788	alg = 't';
1789	store_localhom = 0;
1790	break;
1791	case 'L':
1792	alg = 'L';
1793	break;
1794	case 'Y':
1795	alg = 'Y'; // nadd>0 no toki nomi. moto no hairetsu to atarashii hairetsuno alignmnt -> L;
1796	break;
1797	case 's':
1798	alg = 's';
1799	break;
1800	case 'G':
1801	alg = 'G';
1802	break;
1803	case 'B':
1804	alg = 'B';
1805	break;
1806	case 'T':
1807	alg = 'T';
1808	break;
1809	case 'H':
1810	alg = 'H';
1811	break;
1812	case 'M':
1813	alg = 'M';
1814	break;
1815	case 'R':
1816	alg = 'R';
1817	break;
1818	case 'r':
1819	alg = 'r'; // nadd>0 no toki nomi. moto no hairetsu to atarashii hairetsuno alignmnt -> R, last
1820	break;
1821	case 'N':
1822	alg = 'N';
1823	break;
1824	case 'A':
1825	alg = 'A';
1826	break;
1827	case 'V':
1828	alg = 'V';
1829	break;
1830	case 'F':
1831	use_fft = 1;
1832	break;
1833	case 'v':
1834	tbrweight = 3;
1835	break;
1836	case 'y':
1837	divpairscore = 1;
1838	break;
1839	/* Modified 01/08/27, default: user tree */
1840	case 'J':
1841	tbutree = 0;
1842	break;
1843	/* modification end. */
1844	case 'o':
1845	// foldalignopt = *++argv;
1846	strcat( foldalignopt, " " );
1847	strcat( foldalignopt, *++argv );
1848	fprintf( stderr, "foldalignopt = %s\n", foldalignopt );
1849	--argc;
1850	goto nextoption;
1851	#if 0
1852	case 'z':
1853	fftThreshold = atoi( *++argv );
1854	--argc;
1855	goto nextoption;
1856	case 'w':
1857	fftWinSize = atoi( *++argv );
1858	--argc;
1859	goto nextoption;
1860	case 'Z':
1861	checkC = 1;
1862	break;
1863	#endif
1864	default:
1865	fprintf( stderr, "illegal option %c\n", c );
1866	argc = 0;
1867	break;
1868	}
1869	}
1870	nextoption:
1871	;
1872	}
1873	if( argc == 1 )
1874	{
1875	cut = atof( (*argv) );
1876	argc--;
1877	}
1878	if( argc != 0 )
1879	{
1880	fprintf( stderr, "options: Check source file !\n" );
1881	exit( 1 );
1882	}
1883	if( tbitr == 1 && outgap == 0 )
1884	{
1885	fprintf( stderr, "conflicting options : o, m or u\n" );
1886	exit( 1 );
1887	}
1888	}
1889
1890	int countamino( char *s, int end )
1891	{
1892	int val = 0;
1893	while( end-- )
1894	if( *s++ != '-' ) val++;
1895	return( val );
1896	}
1897
1898	#if enablemultithread
1899	static void athread( void arg ) // alg='R', alg='r' -> tsukawarenai.
1900	{
1901	thread_arg_t targ = (thread_arg_t )arg;
1902	int i, ilim, j, jst;
1903	int off1, off2, dum1, dum2, thereisx;
1904	int intdum;
1905	double bunbo;
1906	float pscore = 0.0; // by D.Mathog
1907	double *effarr1;
1908	double *effarr2;
1909	char mseq1, mseq2, distseq1, distseq2, dumseq1, dumseq2;
1910	char **aseq;
1911
1912	// thread_arg
1913	int thread_no = targ->thread_no;
1914	int njob = targ->njob;
1915	Jobtable *jobpospt = targ->jobpospt;
1916	char **name = targ->name;
1917	char **seq = targ->seq;
1918	char **dseq = targ->dseq;
1919	int *thereisxineachseq = targ->thereisxineachseq;
1920	LocalHom **localhomtable = targ->localhomtable;
1921	double **distancemtx = targ->distancemtx;
1922	double *selfscore = targ->selfscore;
1923	char ***bpp = targ->bpp;
1924	Lastresx **lastresx = targ->lastresx;
1925	int alloclen = targ->alloclen;
1926
1927	// fprintf( stderr, "thread %d start!\n", thread_no );
1928
1929	effarr1 = AllocateDoubleVec( 1 );
1930	effarr2 = AllocateDoubleVec( 1 );
1931	mseq1 = AllocateCharMtx( njob, 0 );
1932	mseq2 = AllocateCharMtx( njob, 0 );
1933	if( alg == 'N' )
1934	{
1935	dumseq1 = AllocateCharMtx( 1, alloclen+10 );
1936	dumseq2 = AllocateCharMtx( 1, alloclen+10 );
1937	}
1938	distseq1 = AllocateCharMtx( 1, 0 );
1939	distseq2 = AllocateCharMtx( 1, 0 );
1940	aseq = AllocateCharMtx( 2, alloclen+10 );
1941
1942	if( alg == 'Y' \|\| alg == 'r' ) ilim = njob - nadd;
1943	else ilim = njob - 1;
1944
1945	while( 1 )
1946	{
1947	pthread_mutex_lock( targ->mutex_counter );
1948	j = jobpospt->j;
1949	i = jobpospt->i;
1950	j++;
1951	if( j == njob )
1952	{
1953	i++;
1954
1955	if( alg == 'Y' \|\| alg == 'r' ) jst = njob - nadd;
1956	else jst = i + 1;
1957	j = jst;
1958
1959	if( i == ilim )
1960	{
1961	// fprintf( stderr, "thread %d end!\n", thread_no );
1962	pthread_mutex_unlock( targ->mutex_counter );
1963
1964	if( commonIP ) FreeIntMtx( commonIP );
1965	commonIP = NULL;
1966	if( commonJP ) FreeIntMtx( commonJP );
1967	commonJP = NULL;
1968	Falign( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL );
1969	G__align11( NULL, NULL, 0, 0, 0 ); // 20130603
1970	G__align11_noalign( NULL, 0, 0, NULL, NULL, 0 );
1971	L__align11( NULL, NULL, 0, NULL, NULL );
1972	L__align11_noalign( NULL, NULL );
1973	genL__align11( NULL, NULL, 0, NULL, NULL );
1974	free( effarr1 );
1975	free( effarr2 );
1976	free( mseq1 );
1977	free( mseq2 );
1978	if( alg == 'N' )
1979	{
1980	FreeCharMtx( dumseq1 );
1981	FreeCharMtx( dumseq2 );
1982	}
1983	free( distseq1 );
1984	free( distseq2 );
1985	FreeCharMtx( aseq );
1986	return( NULL );
1987	}
1988	}
1989	jobpospt->j = j;
1990	jobpospt->i = i;
1991	pthread_mutex_unlock( targ->mutex_counter );
1992
1993
1994	if( j == i+1 \|\| j % 100 == 0 )
1995	{
1996	fprintf( stderr, "% 5d / %d (by thread %3d) \r", i, njob-nadd, thread_no );
1997	// fprintf( stderr, "% 5d - %5d / %d (thread %d)\n", i, j, njob, thread_no );
1998	}
1999
2000
2001	if( strlen( seq[i] ) == 0 \|\| strlen( seq[j] ) == 0 )
2002	{
2003	if( store_dist )
2004	{
2005	if( alg == 'Y' \|\| alg == 'r' ) distancemtx[i][j-(njob-nadd)] = 3.0;
2006	else distancemtx[i][j] = 3.0;
2007	}
2008	if( stdout_dist)
2009	{
2010	pthread_mutex_lock( targ->mutex_stdout );
2011	fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, 3.0 );
2012	pthread_mutex_unlock( targ->mutex_stdout );
2013	}
2014	continue;
2015	}
2016
2017	strcpy( aseq[0], seq[i] );
2018	strcpy( aseq[1], seq[j] );
2019	// clus1 = conjuctionfortbfast( pair, i, aseq, mseq1, effarr1, effarr, indication1 );
2020	// clus2 = conjuctionfortbfast( pair, j, aseq, mseq2, effarr2, effarr, indication2 );
2021	// fprintf( stderr, "Skipping conjuction..\n" );
2022
2023	effarr1[0] = 1.0;
2024	effarr2[0] = 1.0;
2025	mseq1[0] = aseq[0];
2026	mseq2[0] = aseq[1];
2027
2028	thereisx = thereisxineachseq[i] + thereisxineachseq[j];
2029	// strcpy( distseq1[0], dseq[i] ); // nen no tame
2030	// strcpy( distseq2[0], dseq[j] ); // nen no tame
2031	distseq1[0] = dseq[i];
2032	distseq2[0] = dseq[j];
2033
2034	// fprintf( stderr, "mseq1 = %s\n", mseq1[0] );
2035	// fprintf( stderr, "mseq2 = %s\n", mseq2[0] );
2036
2037	#if 0
2038	fprintf( stderr, "group1 = %.66s", indication1 );
2039	fprintf( stderr, "\n" );
2040	fprintf( stderr, "group2 = %.66s", indication2 );
2041	fprintf( stderr, "\n" );
2042	#endif
2043	// for( l=0; l<clus1; l++ ) fprintf( stderr, "## STEP-eff for mseq1-%d %f\n", l, effarr1[l] );
2044
2045	if( use_fft )
2046	{
2047	pscore = Falign( mseq1, mseq2, effarr1, effarr2, 1, 1, alloclen, &intdum, NULL, 0, NULL );
2048	// fprintf( stderr, "pscore (fft) = %f\n", pscore );
2049	off1 = off2 = 0;
2050	}
2051	else
2052	{
2053	switch( alg )
2054	{
2055	case( 'R' ):
2056	if( nadd && njob-nadd <= j && njob-nadd <= i ) // new sequence doushi ha mushi
2057	pscore = 0.0;
2058	else
2059	pscore = (float)lastresx[i][j].score; // all pair
2060	break;
2061	case( 'r' ):
2062	if( nadd == 0 \|\| ( i < njob-nadd && njob-nadd <= j ) )
2063	pscore = (float)lastresx[i][j-(njob-nadd)].score;
2064	else
2065	pscore = 0.0;
2066	break;
2067	case( 'L' ):
2068	if( nadd && njob-nadd <= j && njob-nadd <= i ) // new sequence doushi ha mushi
2069	pscore = 0.0;
2070	else
2071	{
2072	if( store_localhom )
2073	{
2074	pscore = L__align11( mseq1, mseq2, alloclen, &off1, &off2 );
2075	if( thereisx ) pscore = L__align11_noalign( distseq1, distseq2 ); // uwagaki
2076	}
2077	else
2078	pscore = L__align11_noalign( distseq1, distseq2 );
2079	}
2080	// pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, distseq1, distseq2, alloclen ); // CHUUI!!!!!!
2081	break;
2082	case( 'Y' ):
2083	if( nadd == 0 \|\| ( i < njob-nadd && njob-nadd <= j ) ) // new sequence vs exiting sequence nomi keisan
2084	{
2085	if( store_localhom )
2086	{
2087	pscore = L__align11( mseq1, mseq2, alloclen, &off1, &off2 );
2088	if( thereisx ) pscore = L__align11_noalign( distseq1, distseq2 ); // uwagaki
2089	}
2090	else
2091	pscore = L__align11_noalign( distseq1, distseq2 );
2092	}
2093	else
2094	pscore = 0.0;
2095	break;
2096	case( 'A' ):
2097	if( store_localhom )
2098	{
2099	pscore = G__align11( mseq1, mseq2, alloclen, outgap, outgap );
2100	if( thereisx ) pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, distseq1, distseq2, alloclen ); // uwagaki
2101	}
2102	else
2103	pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, distseq1, distseq2, alloclen ); // uwagaki
2104	off1 = off2 = 0;
2105	break;
2106	case( 'N' ):
2107	// pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, alloclen );
2108	pscore = genL__align11( mseq1, mseq2, alloclen, &off1, &off2 );
2109	if( thereisx )
2110	{
2111	strcpy( dumseq1[0], distseq1[0] );
2112	strcpy( dumseq2[0], distseq2[0] );
2113	pscore = genL__align11( dumseq1, dumseq2, alloclen, &dum1, &dum2 ); // uwagaki
2114	}
2115	break;
2116	case( 't' ):
2117	off1 = off2 = 0;
2118	// pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, alloclen );
2119	pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, distseq1, distseq2, alloclen ); // tsuneni distseq shiyou
2120	break;
2121	case( 's' ):
2122	pscore = callmxscarna_giving_bpp( mseq1, mseq2, bpp[i], bpp[j], alloclen, i, j );
2123	off1 = off2 = 0;
2124	break;
2125	case( 'G' ):
2126	pscore = calldafs_giving_bpp( mseq1, mseq2, bpp[i], bpp[j], alloclen, i, j );
2127	off1 = off2 = 0;
2128	break;
2129	#if 0
2130	case( 'a' ):
2131	pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
2132	off1 = off2 = 0;
2133	break;
2134	case( 'K' ):
2135	pscore = genG__align11( mseq1, mseq2, alloclen );
2136	off1 = off2 = 0;
2137	break;
2138	case( 'H' ):
2139	pscore = recallpairfoldalign( mseq1, mseq2, i, j, &off1, &off2, alloclen );
2140	break;
2141	case( 'B' ):
2142	case( 'T' ):
2143	pscore = recalllara( mseq1, mseq2, alloclen );
2144	off1 = off2 = 0;
2145	break;
2146	case( 'M' ):
2147	pscore = MSalign11( mseq1, mseq2, alloclen );
2148	break;
2149	#endif
2150	default:
2151	ErrorExit( "\n\nERROR IN SOURCE FILE\n\n" );
2152	}
2153	}
2154
2155	if( alg == 't' \|\| ( mseq1[0][0] != 0 && mseq2[0][0] != 0 ) ) // 't' no jouken ha iranai to omou. if( ( mseq1[0][0] != 0 && mseq2[0][0] != 0 ) )
2156	{
2157	#if SCOREOUT
2158	fprintf( stderr, "score = %10.2f (%d,%d)\n", pscore, i, j );
2159	#endif
2160	// if( pscore > 0.0 && ( nadd == 0 \|\| ( alg != 'Y' && alg != 'r' ) \|\| ( i < njob-nadd && njob-nadd <= j ) ) ) x-ins-i de seido teika
2161	if( ( nadd == 0 \|\| ( alg != 'Y' && alg != 'r' ) \|\| ( i < njob-nadd && njob-nadd <= j ) ) )
2162	{
2163	if( !store_localhom )
2164	;
2165	else if( alg == 'R' )
2166	putlocalhom_last( mseq1[0], mseq2[0], localhomtable[i]+j, lastresx[i]+j );
2167	else if( alg == 'r' )
2168	putlocalhom_last( mseq1[0], mseq2[0], localhomtable[i]+j-(njob-nadd), lastresx[i]+j-(njob-nadd) );// ?????
2169	else if( alg == 'H' )
2170	putlocalhom_ext( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
2171	else if( alg == 'Y' )
2172	putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j-(njob-nadd), off1, off2, (int)pscore, strlen( mseq1[0] ) );
2173	else if( alg != 'S' && alg != 'V' )
2174	putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
2175	}
2176
2177	if( (bunbo=MIN( selfscore[i], selfscore[j] )) == 0.0 )
2178	pscore = 2.0;
2179	else if( bunbo < pscore ) // mondai ari
2180	pscore = 0.0;
2181	else
2182	pscore = ( 1.0 - pscore / bunbo ) * 2.0;
2183	}
2184	else
2185	{
2186	pscore = 2.0;
2187	}
2188
2189	#if 1 // mutex
2190	if( stdout_align )
2191	{
2192	pthread_mutex_lock( targ->mutex_stdout );
2193	if( alg != 't' )
2194	{
2195	fprintf( stdout, "sequence %d - sequence %d, pairwise distance = %10.5f\n", i+1, j+1, pscore );
2196	fprintf( stdout, ">%s\n", name[i] );
2197	write1seq( stdout, mseq1[0] );
2198	fprintf( stdout, ">%s\n", name[j] );
2199	write1seq( stdout, mseq2[0] );
2200	fprintf( stdout, "\n" );
2201	}
2202	pthread_mutex_unlock( targ->mutex_stdout );
2203	}
2204	if( stdout_dist )
2205	{
2206	pthread_mutex_lock( targ->mutex_stdout );
2207	if( j == i+1 ) fprintf( stdout, "%d %d d=%.3f\n", i+1, i+1, 0.0 );
2208	fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, pscore );
2209	pthread_mutex_unlock( targ->mutex_stdout );
2210	}
2211	#endif // mutex
2212	if( store_dist )
2213	{
2214	if( alg == 'Y' \|\| alg == 'r' ) distancemtx[i][j-(njob-nadd)] = pscore;
2215	else distancemtx[i][j] = pscore;
2216	}
2217	}
2218	}
2219	#endif
2220
2221	static void pairalign( char *name, int nlen, char seq, char aseq, char *dseq, int thereisxineachseq, char mseq1, char mseq2, int alloclen, Lastresx **lastresx )
2222	{
2223	int i, j, ilim, jst, jj;
2224	int off1, off2, dum1, dum2, thereisx;
2225	float pscore = 0.0; // by D.Mathog
2226	FILE hat2p, hat3p;
2227	double **distancemtx;
2228	double *selfscore;
2229	double *effarr1;
2230	double *effarr2;
2231	char *pt;
2232	char *hat2file = "hat2";
2233	LocalHom *localhomtable = NULL, tmpptr;
2234	int intdum;
2235	double bunbo;
2236	char ***bpp = NULL; // mxscarna no toki dake
2237	char distseq1, distseq2;
2238	char dumseq1, dumseq2;
2239
2240	if( store_localhom )
2241	{
2242	if( alg == 'Y' \|\| alg == 'r' )
2243	{
2244	ilim = njob - nadd;
2245	jst = nadd;
2246	}
2247	else
2248	{
2249	ilim = njob;
2250	jst = njob;
2251	}
2252	localhomtable = (LocalHom *)calloc( ilim, sizeof( LocalHom ) );
2253	for( i=0; i<ilim; i++)
2254	{
2255	localhomtable[i] = (LocalHom *)calloc( jst, sizeof( LocalHom ) );
2256	for( j=0; j<jst; j++)
2257	{
2258	localhomtable[i][j].start1 = -1;
2259	localhomtable[i][j].end1 = -1;
2260	localhomtable[i][j].start2 = -1;
2261	localhomtable[i][j].end2 = -1;
2262	localhomtable[i][j].opt = -1.0;
2263	localhomtable[i][j].next = NULL;
2264	localhomtable[i][j].nokori = 0;
2265	}
2266	}
2267	}
2268
2269	if( store_dist )
2270	{
2271	if( alg == 'Y' \|\| alg == 'r' )
2272	distancemtx = AllocateDoubleMtx( njob, nadd );
2273	else
2274	distancemtx = AllocateDoubleMtx( njob, njob );
2275	}
2276	else distancemtx = NULL;
2277
2278	if( alg == 'N' )
2279	{
2280	dumseq1 = AllocateCharMtx( 1, alloclen+10 );
2281	dumseq2 = AllocateCharMtx( 1, alloclen+10 );
2282	}
2283	distseq1 = AllocateCharMtx( 1, 0 ); // muda
2284	distseq2 = AllocateCharMtx( 1, 0 ); // muda
2285
2286	selfscore = AllocateDoubleVec( njob );
2287	effarr1 = AllocateDoubleVec( njob );
2288	effarr2 = AllocateDoubleVec( njob );
2289
2290	#if 0
2291	fprintf( stderr, "##### fftwinsize = %d, fftthreshold = %d\n", fftWinSize, fftThreshold );
2292	#endif
2293
2294	#if 0
2295	for( i=0; i<njob; i++ )
2296	fprintf( stderr, "TBFAST effarr[%d] = %f\n", i, effarr[i] );
2297	#endif
2298
2299
2300	// writePre( njob, name, nlen, aseq, 0 );
2301
2302	if( alg == 'R' )
2303	{
2304	fprintf( stderr, "Calling last (http://last.cbrc.jp/)\n" );
2305	if( lastonce )
2306	calllast_once( njob, seq, njob, seq, lastresx );
2307	else
2308	calllast_fast( njob, seq, njob, seq, lastresx );
2309	fprintf( stderr, "done.\n" );
2310	// nthread = 0; // igo multithread nashi
2311	}
2312	if( alg == 'r' )
2313	{
2314	fprintf( stderr, "Calling last (http://last.cbrc.jp/)\n" );
2315	fprintf( stderr, "nadd=%d\n", nadd );
2316	#if 1 // last_fast ha, lastdb ga muda
2317	if( lastonce )
2318	calllast_once( njob-nadd, seq, nadd, seq+njob-nadd, lastresx );
2319	else
2320	calllast_fast( njob-nadd, seq, nadd, seq+njob-nadd, lastresx );
2321	#else
2322	calllast_once( njob-nadd, seq, nadd, seq+njob-nadd, lastresx );
2323	#endif
2324
2325	fprintf( stderr, "nadd=%d\n", nadd );
2326	fprintf( stderr, "done.\n" );
2327	// nthread = 0; // igo multithread nashi
2328	}
2329
2330	if( alg == 'H' )
2331	{
2332	fprintf( stderr, "Calling FOLDALIGN with option '%s'\n", foldalignopt );
2333	callfoldalign( njob, seq );
2334	fprintf( stderr, "done.\n" );
2335	}
2336	if( alg == 'B' )
2337	{
2338	fprintf( stderr, "Running LARA (Bauer et al. http://www.planet-lisa.net/)\n" );
2339	calllara( njob, seq, "" );
2340	fprintf( stderr, "done.\n" );
2341	}
2342	if( alg == 'T' )
2343	{
2344	fprintf( stderr, "Running SLARA (Bauer et al. http://www.planet-lisa.net/)\n" );
2345	calllara( njob, seq, "-s" );
2346	fprintf( stderr, "done.\n" );
2347	}
2348	if( alg == 's' )
2349	{
2350	fprintf( stderr, "Preparing bpp\n" );
2351	// bpp = AllocateCharCub( njob, nlenmax, 0 );
2352	bpp = calloc( njob, sizeof( char ** ) );
2353	preparebpp( njob, bpp );
2354	fprintf( stderr, "done.\n" );
2355	fprintf( stderr, "Running MXSCARNA (Tabei et al. http://www.ncrna.org/software/mxscarna)\n" );
2356	}
2357	if( alg == 'G' )
2358	{
2359	fprintf( stderr, "Preparing bpp\n" );
2360	// bpp = AllocateCharCub( njob, nlenmax, 0 );
2361	bpp = calloc( njob, sizeof( char ** ) );
2362	preparebpp( njob, bpp );
2363	fprintf( stderr, "done.\n" );
2364	fprintf( stderr, "Running DAFS (Sato et al. http://www.ncrna.org/)\n" );
2365	}
2366
2367	for( i=0; i<njob; i++ )
2368	{
2369	pscore = 0.0;
2370	for( pt=seq[i]; *pt; pt++ )
2371	pscore += amino_dis[(int)pt][(int)pt];
2372	selfscore[i] = pscore;
2373	// fprintf( stderr, "selfscore[%d] = %f\n", i, selfscore[i] );
2374	}
2375
2376	#if enablemultithread
2377	if( nthread > 0 ) // alg=='r' \|\| alg=='R' -> nthread:=0 (sukoshi ue)
2378	{
2379	Jobtable jobpos;
2380	pthread_t *handle;
2381	pthread_mutex_t mutex_counter;
2382	pthread_mutex_t mutex_stdout;
2383	thread_arg_t *targ;
2384
2385	if( alg == 'Y' \|\| alg == 'r' ) jobpos.j = njob - nadd - 1;
2386	else jobpos.j = 0;
2387	jobpos.i = 0;
2388
2389	targ = calloc( nthread, sizeof( thread_arg_t ) );
2390	handle = calloc( nthread, sizeof( pthread_t ) );
2391	pthread_mutex_init( &mutex_counter, NULL );
2392	pthread_mutex_init( &mutex_stdout, NULL );
2393
2394	for( i=0; i<nthread; i++ )
2395	{
2396	targ[i].thread_no = i;
2397	targ[i].njob = njob;
2398	targ[i].jobpospt = &jobpos;
2399	targ[i].name = name;
2400	targ[i].seq = seq;
2401	targ[i].dseq = dseq;
2402	targ[i].thereisxineachseq = thereisxineachseq;
2403	targ[i].localhomtable = localhomtable;
2404	targ[i].distancemtx = distancemtx;
2405	targ[i].selfscore = selfscore;
2406	targ[i].bpp = bpp;
2407	targ[i].lastresx = lastresx;
2408	targ[i].alloclen = alloclen;
2409	targ[i].mutex_counter = &mutex_counter;
2410	targ[i].mutex_stdout = &mutex_stdout;
2411
2412	// athread( (void *)targ );
2413	pthread_create( handle+i, NULL, athread, (void *)(targ+i) );
2414	// pthread_create( handle+i, NULL, bthread, (void *)(targ+i) );
2415	}
2416
2417
2418	for( i=0; i<nthread; i++ )
2419	{
2420	pthread_join( handle[i], NULL );
2421	}
2422	pthread_mutex_destroy( &mutex_counter );
2423	pthread_mutex_destroy( &mutex_stdout );
2424	free( handle );
2425	free( targ );
2426	}
2427	else
2428	#endif
2429	{
2430	if( alg == 'Y' \|\| alg == 'r' ) ilim = njob - nadd;
2431	else ilim = njob - 1;
2432	for( i=0; i<ilim; i++ )
2433	{
2434	if( stdout_dist) fprintf( stdout, "%d %d d=%.3f\n", i+1, i+1, 0.0 );
2435	fprintf( stderr, "% 5d / %d\r", i, njob-nadd );
2436	fflush( stderr );
2437
2438	if( alg == 'Y' \|\| alg == 'r' ) jst = njob - nadd;
2439	else jst = i + 1;
2440	for( j=jst; j<njob; j++ )
2441	{
2442
2443	if( strlen( seq[i] ) == 0 \|\| strlen( seq[j] ) == 0 )
2444	{
2445	if( store_dist )
2446	{
2447	if( alg == 'Y' \|\| alg == 'r' ) distancemtx[i][j-(njob-nadd)] = 3.0;
2448	else distancemtx[i][j] = 3.0;
2449	}
2450	if( stdout_dist) fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, 3.0 );
2451	continue;
2452	}
2453
2454	strcpy( aseq[0], seq[i] );
2455	strcpy( aseq[1], seq[j] );
2456	// clus1 = conjuctionfortbfast( pair, i, aseq, mseq1, effarr1, effarr, indication1 );
2457	// clus2 = conjuctionfortbfast( pair, j, aseq, mseq2, effarr2, effarr, indication2 );
2458	// fprintf( stderr, "Skipping conjuction..\n" );
2459
2460	effarr1[0] = 1.0;
2461	effarr2[0] = 1.0;
2462	mseq1[0] = aseq[0];
2463	mseq2[0] = aseq[1];
2464
2465	thereisx = thereisxineachseq[i] + thereisxineachseq[j];
2466	// strcpy( distseq1[0], dseq[i] ); // nen no tame
2467	// strcpy( distseq2[0], dseq[j] ); // nen no tame
2468	distseq1[0] = dseq[i];
2469	distseq2[0] = dseq[j];
2470
2471	// fprintf( stderr, "mseq1 = %s\n", mseq1[0] );
2472	// fprintf( stderr, "mseq2 = %s\n", mseq2[0] );
2473
2474	#if 0
2475	fprintf( stderr, "group1 = %.66s", indication1 );
2476	fprintf( stderr, "\n" );
2477	fprintf( stderr, "group2 = %.66s", indication2 );
2478	fprintf( stderr, "\n" );
2479	#endif
2480	// for( l=0; l<clus1; l++ ) fprintf( stderr, "## STEP-eff for mseq1-%d %f\n", l, effarr1[l] );
2481
2482	if( use_fft )
2483	{
2484	pscore = Falign( mseq1, mseq2, effarr1, effarr2, 1, 1, alloclen, &intdum, NULL, 0, NULL );
2485	// fprintf( stderr, "pscore (fft) = %f\n", pscore );
2486	off1 = off2 = 0;
2487	}
2488	else
2489	{
2490	switch( alg )
2491	{
2492	case( 't' ):
2493	// pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, alloclen );
2494	pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, distseq1, distseq2, alloclen ); // tsuneni distseq shiyou
2495	off1 = off2 = 0;
2496	break;
2497	case( 'A' ):
2498	if( store_localhom )
2499	{
2500	pscore = G__align11( mseq1, mseq2, alloclen, outgap, outgap );
2501	if( thereisx ) pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, distseq1, distseq2, alloclen ); // uwagaki
2502	}
2503	else
2504	pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, distseq1, distseq2, alloclen ); // uwagaki
2505	off1 = off2 = 0;
2506	break;
2507	case( 'N' ):
2508	// pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, alloclen );
2509	pscore = genL__align11( mseq1, mseq2, alloclen, &off1, &off2 );
2510	if( thereisx )
2511	{
2512	strcpy( dumseq1[0], distseq1[0] );
2513	strcpy( dumseq2[0], distseq2[0] );
2514	pscore = genL__align11( dumseq1, dumseq2, alloclen, &dum1, &dum2 ); // uwagaki
2515	}
2516	break;
2517	case( 'R' ):
2518	if( nadd && njob-nadd <= j && njob-nadd <= i ) // new sequence doushi ha mushi
2519	pscore = 0.0;
2520	else
2521	pscore = (float)lastresx[i][j].score; // all pair
2522	break;
2523	case( 'r' ):
2524	if( nadd == 0 \|\| ( i < njob-nadd && njob-nadd <= j ) )
2525	pscore = (float)lastresx[i][j-(njob-nadd)].score;
2526	else
2527	pscore = 0.0;
2528	break;
2529	case( 'L' ):
2530	if( nadd && njob-nadd <= j && njob-nadd <= i ) // new sequence doushi ha mushi
2531	pscore = 0.0;
2532	else
2533	{
2534	if( store_localhom )
2535	{
2536	pscore = L__align11( mseq1, mseq2, alloclen, &off1, &off2 ); // all pair
2537	if( thereisx ) pscore = L__align11_noalign( distseq1, distseq2 ); // all pair
2538	}
2539	else
2540	pscore = L__align11_noalign( distseq1, distseq2 ); // all pair
2541	}
2542	// pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, distseq1, distseq2, alloclen ); // CHUUI!!!!!!
2543	break;
2544	case( 'Y' ):
2545	if( nadd == 0 \|\| ( i < njob-nadd && njob-nadd <= j ) ) // new sequence vs exiting sequence nomi keisan
2546	{
2547	if( store_localhom )
2548	{
2549	pscore = L__align11( mseq1, mseq2, alloclen, &off1, &off2 );
2550	if( thereisx ) pscore = L__align11_noalign( distseq1, distseq2 ); // uwagaki
2551	}
2552	else
2553	pscore = L__align11_noalign( distseq1, distseq2 );
2554	}
2555	else
2556	pscore = 0.0;
2557	break;
2558	case( 'a' ):
2559	pscore = Aalign( mseq1, mseq2, effarr1, effarr2, 1, 1, alloclen );
2560	off1 = off2 = 0;
2561	break;
2562	#if 0
2563	case( 'K' ):
2564	pscore = genG__align11( mseq1, mseq2, alloclen );
2565	off1 = off2 = 0;
2566	break;
2567	#endif
2568	case( 'H' ):
2569	pscore = recallpairfoldalign( mseq1, mseq2, i, j, &off1, &off2, alloclen );
2570	break;
2571	case( 'B' ):
2572	case( 'T' ):
2573	pscore = recalllara( mseq1, mseq2, alloclen );
2574	off1 = off2 = 0;
2575	break;
2576	case( 's' ):
2577	pscore = callmxscarna_giving_bpp( mseq1, mseq2, bpp[i], bpp[j], alloclen, i, j );
2578	off1 = off2 = 0;
2579	break;
2580	case( 'G' ):
2581	pscore = calldafs_giving_bpp( mseq1, mseq2, bpp[i], bpp[j], alloclen, i, j );
2582	off1 = off2 = 0;
2583	break;
2584	case( 'M' ):
2585	pscore = MSalign11( mseq1, mseq2, alloclen );
2586	break;
2587	default:
2588	ErrorExit( "ERROR IN SOURCE FILE" );
2589	}
2590	}
2591
2592	if( alg == 't' \|\| ( mseq1[0][0] != 0 && mseq2[0][0] != 0 ) ) // 't' no jouken ha iranai to omou. if( ( mseq1[0][0] != 0 && mseq2[0][0] != 0 ) )
2593	{
2594	#if SCOREOUT
2595	fprintf( stderr, "score = %10.2f (%d,%d)\n", pscore, i, j );
2596	#endif
2597	// if( pscore > 0.0 && ( nadd == 0 \|\| ( alg != 'Y' && alg != 'r' ) \|\| ( i < njob-nadd && njob-nadd <= j ) ) ) // x-ins-i de seido teika
2598	if( ( nadd == 0 \|\| ( alg != 'Y' && alg != 'r' ) \|\| ( i < njob-nadd && njob-nadd <= j ) ) )
2599	{
2600	if( !store_localhom )
2601	;
2602	else if( alg == 'R' )
2603	putlocalhom_last( mseq1[0], mseq2[0], localhomtable[i]+j, lastresx[i]+j );
2604	else if( alg == 'r' )
2605	putlocalhom_last( mseq1[0], mseq2[0], localhomtable[i]+j-(njob-nadd), lastresx[i]+j-(njob-nadd) );// ?????
2606	else if( alg == 'H' )
2607	putlocalhom_ext( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
2608	else if( alg == 'Y' )
2609	putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j-(njob-nadd), off1, off2, (int)pscore, strlen( mseq1[0] ) );
2610	else if( alg != 'S' && alg != 'V' )
2611	putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
2612	}
2613
2614	if( (bunbo=MIN( selfscore[i], selfscore[j] )) == 0.0 )
2615	pscore = 2.0;
2616	else if( bunbo < pscore ) // mondai ari
2617	pscore = 0.0;
2618	else
2619	pscore = ( 1.0 - pscore / bunbo ) * 2.0;
2620	}
2621	else
2622	{
2623	pscore = 2.0;
2624	}
2625
2626	if( stdout_align )
2627	{
2628	if( alg != 't' )
2629	{
2630	fprintf( stdout, "sequence %d - sequence %d, pairwise distance = %10.5f\n", i+1, j+1, pscore );
2631	fprintf( stdout, ">%s\n", name[i] );
2632	write1seq( stdout, mseq1[0] );
2633	fprintf( stdout, ">%s\n", name[j] );
2634	write1seq( stdout, mseq2[0] );
2635	fprintf( stdout, "\n" );
2636	}
2637	}
2638	if( stdout_dist ) fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, pscore );
2639	if( store_dist)
2640	{
2641	if( alg == 'Y' \|\| alg == 'r' ) distancemtx[i][j-(njob-nadd)] = pscore;
2642	else distancemtx[i][j] = pscore;
2643	}
2644	}
2645	}
2646	}
2647
2648
2649	if( store_dist )
2650	{
2651	hat2p = fopen( hat2file, "w" );
2652	if( !hat2p ) ErrorExit( "Cannot open hat2." );
2653	if( alg == 'Y' \|\| alg == 'r' )
2654	WriteHat2_part_pointer( hat2p, njob, nadd, name, distancemtx );
2655	else
2656	WriteHat2_pointer( hat2p, njob, name, distancemtx );
2657	fclose( hat2p );
2658	}
2659
2660	hat3p = fopen( "hat3", "w" );
2661	if( !hat3p ) ErrorExit( "Cannot open hat3." );
2662	if( store_localhom )
2663	{
2664	fprintf( stderr, "\n\n##### writing hat3\n" );
2665	if( alg == 'Y' \|\| alg == 'r' )
2666	ilim = njob-nadd;
2667	else
2668	ilim = njob-1;
2669	for( i=0; i<ilim; i++ )
2670	{
2671	if( alg == 'Y' \|\| alg == 'r' )
2672	{
2673	jst = njob-nadd;
2674	jj = 0;
2675	}
2676	else
2677	{
2678	jst = i+1;
2679	jj = i+1;
2680	}
2681	for( j=jst; j<njob; j++, jj++ )
2682	{
2683	for( tmpptr=localhomtable[i]+jj; tmpptr; tmpptr=tmpptr->next )
2684	{
2685	// fprintf( stderr, "j=%d, jj=%d\n", j, jj );
2686	if( tmpptr->opt == -1.0 ) continue;
2687	// tmptmptmptmptmp
2688	// if( alg == 'B' \|\| alg == 'T' )
2689	// fprintf( hat3p, "%d %d %d %7.5f %d %d %d %d %p\n", i, j, tmpptr->overlapaa, 1.0, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2, (void *)tmpptr->next );
2690	// else
2691	fprintf( hat3p, "%d %d %d %7.5f %d %d %d %d h\n", i, j, tmpptr->overlapaa, tmpptr->opt, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2 );
2692	// fprintf( hat3p, "%d %d %d %7.5f %d %d %d %d h\n", i, j, tmpptr->overlapaa, tmpptr->opt, tmpptr->start1, tmpptr->end1, tmpptr->start2+1, tmpptr->end2+1 ); // zettai dame!!!!
2693	}
2694	}
2695	}
2696	#if DEBUG
2697	fprintf( stderr, "calling FreeLocalHomTable\n" );
2698	#endif
2699	if( alg == 'Y' \|\| alg == 'r' )
2700	FreeLocalHomTable_part( localhomtable, (njob-nadd), nadd );
2701	else
2702	FreeLocalHomTable( localhomtable, njob );
2703	#if DEBUG
2704	fprintf( stderr, "done. FreeLocalHomTable\n" );
2705	#endif
2706	}
2707	fclose( hat3p );
2708
2709	if( alg == 's' )
2710	{
2711	char **ptpt;
2712	for( i=0; i<njob; i++ )
2713	{
2714	ptpt = bpp[i];
2715	while( 1 )
2716	{
2717	if( ptpt ) free( ptpt );
2718	else break;
2719	ptpt++;
2720	}
2721	free( bpp[i] );
2722	}
2723	free( bpp );
2724	}
2725	free( selfscore );
2726	free( effarr1 );
2727	free( effarr2 );
2728	if( alg == 'N' )
2729	{
2730	FreeCharMtx( dumseq1 );
2731	FreeCharMtx( dumseq2 );
2732	}
2733	free( distseq1 );
2734	free( distseq2 );
2735	if( store_dist ) FreeDoubleMtx( distancemtx );
2736	}
2737
2738
2739
2740	int main( int argc, char *argv[] )
2741	{
2742	int nlen, thereisxineachseq;
2743	char name, seq;
2744	char mseq1, mseq2;
2745	char **aseq;
2746	char **bseq;
2747	char **dseq;
2748	int i, j, k;
2749	FILE *infp;
2750	char c;
2751	int alloclen;
2752	Lastresx **lastresx;
2753
2754	arguments( argc, argv );
2755	#ifndef enablemultithread
2756	nthread = 0;
2757	#endif
2758	fprintf( stderr, "lastonce = %d\n", lastonce );
2759
2760	if( inputfile )
2761	{
2762	infp = fopen( inputfile, "r" );
2763	if( !infp )
2764	{
2765	fprintf( stderr, "Cannot open %s\n", inputfile );
2766	exit( 1 );
2767	}
2768	}
2769	else
2770	infp = stdin;
2771
2772	getnumlen( infp );
2773	rewind( infp );
2774
2775	if( njob < 2 )
2776	{
2777	fprintf( stderr, "At least 2 sequences should be input!\n"
2778	"Only %d sequence found.\n", njob );
2779	exit( 1 );
2780	}
2781	if( njob > M )
2782	{
2783	fprintf( stderr, "The number of sequences must be < %d\n", M );
2784	fprintf( stderr, "Please try the splittbfast program for such large data.\n" );
2785	exit( 1 );
2786	}
2787
2788	if( ( alg == 'r' \|\| alg == 'R' ) && dorp == 'p' )
2789	{
2790	fprintf( stderr, "Not yet supported\n" );
2791	exit( 1 );
2792	}
2793
2794	alloclen = nlenmax*2;
2795	seq = AllocateCharMtx( njob, alloclen+10 );
2796	aseq = AllocateCharMtx( 2, alloclen+10 );
2797	bseq = AllocateCharMtx( njob, alloclen+10 );
2798	dseq = AllocateCharMtx( njob, alloclen+10 );
2799	mseq1 = AllocateCharMtx( njob, 0 );
2800	mseq2 = AllocateCharMtx( njob, 0 );
2801	name = AllocateCharMtx( njob, B );
2802	nlen = AllocateIntVec( njob );
2803	thereisxineachseq = AllocateIntVec( njob );
2804
2805	if( alg == 'R' )
2806	{
2807	lastresx = calloc( njob+1, sizeof( Lastresx * ) );
2808	for( i=0; i<njob; i++ )
2809	{
2810	lastresx[i] = calloc( njob+1, sizeof( Lastresx ) ); // muda
2811	for( j=0; j<njob; j++ )
2812	{
2813	lastresx[i][j].score = 0;
2814	lastresx[i][j].naln = 0;
2815	lastresx[i][j].aln = NULL;
2816	}
2817	lastresx[i][njob].naln = -1;
2818	}
2819	lastresx[njob] = NULL;
2820	}
2821	else if( alg == 'r' )
2822	{
2823	// fprintf( stderr, "Allocating lastresx (%d), njob=%d, nadd=%d\n", njob-nadd+1, njob, nadd );
2824	lastresx = calloc( njob-nadd+1, sizeof( Lastresx * ) );
2825	for( i=0; i<njob-nadd; i++ )
2826	{
2827	// fprintf( stderr, "Allocating lastresx[%d]\n", i );
2828	lastresx[i] = calloc( nadd+1, sizeof( Lastresx ) );
2829	for( j=0; j<nadd; j++ )
2830	{
2831	// fprintf( stderr, "Initializing lastresx[%d][%d]\n", i, j );
2832	lastresx[i][j].score = 0;
2833	lastresx[i][j].naln = 0;
2834	lastresx[i][j].aln = NULL;
2835	}
2836	lastresx[i][nadd].naln = -1;
2837	}
2838	lastresx[njob-nadd] = NULL;
2839	}
2840	else
2841	lastresx = NULL;
2842
2843	#if 0
2844	Read( name, nlen, seq );
2845	#else
2846	readData_pointer( infp, name, nlen, seq );
2847	#endif
2848	fclose( infp );
2849
2850	constants( njob, seq );
2851
2852	#if 0
2853	fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
2854	#endif
2855
2856	initSignalSM();
2857
2858	initFiles();
2859
2860	// WriteOptions( trap_g );
2861
2862	c = seqcheck( seq );
2863	if( c )
2864	{
2865	fprintf( stderr, "Illegal character %c\n", c );
2866	exit( 1 );
2867	}
2868
2869	// writePre( njob, name, nlen, seq, 0 );
2870
2871
2872	for( i=0; i<njob; i++ )
2873	{
2874	gappick0( bseq[i], seq[i] );
2875	thereisxineachseq[i] = removex( dseq[i], bseq[i] );
2876	}
2877
2878	pairalign( name, nlen, bseq, aseq, dseq, thereisxineachseq, mseq1, mseq2, alloclen, lastresx );
2879
2880	fprintf( trap_g, "done.\n" );
2881	#if DEBUG
2882	fprintf( stderr, "closing trap_g\n" );
2883	#endif
2884	fclose( trap_g );
2885
2886	// writePre( njob, name, nlen, aseq, !contin );
2887	#if 0
2888	writeData( stdout, njob, name, nlen, aseq );
2889	#endif
2890	#if IODEBUG
2891	fprintf( stderr, "OSHIMAI\n" );
2892	#endif
2893	SHOWVERSION;
2894
2895	if( stdout_dist && nthread > 1 )
2896	{
2897	fprintf( stderr, "\nThe order of distances is not identical to that in the input file, because of the parallel calculation. Reorder them by yourself, using sort -n -k 2 \| sort -n -k 1 -s\n" );
2898	}
2899	if( stdout_align && nthread > 1 )
2900	{
2901	fprintf( stderr, "\nThe order of pairwise alignments is not identical to that in the input file, because of the parallel calculation. Reorder them by yourself.\n" );
2902	}
2903
2904	#if 1
2905	if( lastresx )
2906	{
2907	for( i=0; lastresx[i]; i++ )
2908	{
2909	for( j=0; lastresx[i][j].naln!=-1; j++ )
2910	{
2911	for( k=0; k<lastresx[i][j].naln; k++ )
2912	{
2913	free( lastresx[i][j].aln[k].reg1 );
2914	free( lastresx[i][j].aln[k].reg2 );
2915	}
2916	free( lastresx[i][j].aln );
2917	}
2918	free( lastresx[i] );
2919	}
2920	free( lastresx );
2921	}
2922	#endif
2923	FreeCharMtx( seq );
2924	FreeCharMtx( aseq );
2925	FreeCharMtx( bseq );
2926	FreeCharMtx( dseq );
2927	FreeCharMtx( name );
2928	free( mseq1 );
2929	free( mseq2 );
2930	free( nlen );
2931	free( thereisxineachseq );
2932
2933	return( 0 );
2934	}

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

Context Navigation

source: branches/items/GDE/MAFFT/mafft-7.055-with-extensions/core/pairlocalalign.c

Download in other formats: