Context Navigation

makedirectionlist.c.last

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: 19.2 KB

Line
1	#include "mltaln.h"
2
3	#define DEBUG 0
4	#define IODEBUG 0
5	#define SCOREOUT 0
6
7	#define END_OF_VEC -1
8
9	int nadd;
10	float thresholdtorev;
11	int dodp;
12	int addfragment;
13	static char *whereispairalign;
14	static int laste;
15	static int lastm;
16
17	typedef struct _thread_arg
18	{
19	int iend;
20	char **seq;
21	char *tmpseq;
22	int *res;
23	int **spointt;
24	short *table1;
25	int iq;
26	#ifdef enablemultithread
27	int *jshare;
28	int thread_no;
29	pthread_mutex_t *mutex_counter;
30	#endif
31	} thread_arg_t;
32
33
34
35	void arguments( int argc, char *argv[] )
36	{
37	int c;
38
39	nthread = 1;
40	inputfile = NULL;
41	nadd = 0;
42	dodp = 0;
43	alg = 'a';
44	alg = 'm';
45	dorp = NOTSPECIFIED;
46	fmodel = 0;
47	// ppenalty = (int)( -2.0 * 1000 - 0.5 );
48	// ppenalty_ex = (int)( -0.1 * 1000 - 0.5 );
49	// poffset = (int)( 0.1 * 1000 - 0.5 );
50	ppenalty = NOTSPECIFIED;
51	ppenalty_ex = NOTSPECIFIED;
52	poffset = NOTSPECIFIED;
53	kimuraR = 2;
54	pamN = 200;
55	thresholdtorev = 0.1;
56	addfragment = 0;
57	laste = 5000;
58	lastm = 3;
59
60
61	while( --argc > 0 && (*++argv)[0] == '-' )
62	{
63	while ( (c = *++argv[0]) )
64	{
65	switch( c )
66	{
67	case 'i':
68	inputfile = *++argv;
69	fprintf( stderr, "inputfile = %s\n", inputfile );
70	--argc;
71	goto nextoption;
72	case 'I':
73	nadd = atoi( *++argv );
74	fprintf( stderr, "nadd = %d\n", nadd );
75	--argc;
76	goto nextoption;
77	case 'C':
78	nthread = atoi( *++argv );
79	fprintf( stderr, "nthread = %d\n", nthread );
80	--argc;
81	goto nextoption;
82	case 'f':
83	ppenalty = (int)( atof( ++argv ) 1000 - 0.5 );
84	// fprintf( stderr, "ppenalty = %d\n", ppenalty );
85	--argc;
86	goto nextoption;
87	case 'g':
88	ppenalty_ex = (int)( atof( ++argv ) 1000 - 0.5 );
89	fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
90	--argc;
91	goto nextoption;
92	case 'h':
93	poffset = (int)( atof( ++argv ) 1000 - 0.5 );
94	// fprintf( stderr, "poffset = %d\n", poffset );
95	--argc;
96	goto nextoption;
97	case 'k':
98	kimuraR = atoi( *++argv );
99	fprintf( stderr, "kappa = %d\n", kimuraR );
100	--argc;
101	goto nextoption;
102	case 'j':
103	pamN = atoi( *++argv );
104	scoremtx = 0;
105	TMorJTT = JTT;
106	fprintf( stderr, "jtt/kimura %d\n", pamN );
107	--argc;
108	goto nextoption;
109	case 't':
110	thresholdtorev = atof( *++argv );
111	fprintf( stderr, "thresholdtorev = %f\n", thresholdtorev );
112	--argc;
113	goto nextoption;
114	case 'D':
115	whereispairalign = *++argv;
116	fprintf( stderr, "whereispairalign = %s\n", whereispairalign );
117	--argc;
118	goto nextoption;
119	case 'w':
120	lastm = atoi( *++argv );
121	fprintf( stderr, "lastm = %d\n", lastm );
122	--argc;
123	goto nextoption;
124	case 'e':
125	laste = atoi( *++argv );
126	fprintf( stderr, "laste = %d\n", laste );
127	--argc;
128	goto nextoption;
129	case 'd':
130	dodp = 1;
131	break;
132	case 'l':
133	dodp = 2;
134	break;
135	case 'F':
136	addfragment = 1;
137	break;
138	#if 1
139	case 'a':
140	fmodel = 1;
141	break;
142	#endif
143	case 'S':
144	alg = 'S';
145	break;
146	case 'M':
147	alg = 'M';
148	break;
149	case 'm':
150	alg = 'm';
151	break;
152	case 'G':
153	alg = 'G';
154	break;
155	case 'P':
156	dorp = 'p';
157	break;
158	default:
159	fprintf( stderr, "illegal option %c\n", c );
160	argc = 0;
161	break;
162	}
163	}
164	nextoption:
165	;
166	}
167	if( argc == 1 )
168	{
169	cut = atof( (*argv) );
170	argc--;
171	}
172	if( argc != 0 )
173	{
174	fprintf( stderr, "options: Check source file !\n" );
175	exit( 1 );
176	}
177	if( tbitr == 1 && outgap == 0 )
178	{
179	fprintf( stderr, "conflicting options : o, m or u\n" );
180	exit( 1 );
181	}
182	}
183
184
185
186
187	static int maxl;
188	static int tsize;
189
190	void seq_grp_nuc( int grp, char seq )
191	{
192	int tmp;
193	int *grpbk = grp;
194	while( *seq )
195	{
196	tmp = amino_grp[(int)*seq++];
197	if( tmp < 4 )
198	*grp++ = tmp;
199	else
200	// fprintf( stderr, "WARNING : Unknown character %c\r", *(seq-1) );
201	;
202	}
203	*grp = END_OF_VEC;
204	if( grp - grpbk < 6 )
205	{
206	// fprintf( stderr, "\n\nWARNING: Too short.\nPlease also consider use mafft-ginsi, mafft-linsi or mafft-ginsi.\n\n\n" );
207	// exit( 1 );
208	*grpbk = -1;
209	}
210	}
211
212	void seq_grp( int grp, char seq )
213	{
214	int tmp;
215	int *grpbk = grp;
216	while( *seq )
217	{
218	tmp = amino_grp[(int)*seq++];
219	if( tmp < 6 )
220	*grp++ = tmp;
221	else
222	// fprintf( stderr, "WARNING : Unknown character %c\r", *(seq-1) );
223	;
224	}
225	*grp = END_OF_VEC;
226	if( grp - grpbk < 6 )
227	{
228	// fprintf( stderr, "\n\nWARNING: Too short.\nPlease also consider use mafft-ginsi, mafft-linsi or mafft-ginsi.\n\n\n" );
229	// exit( 1 );
230	*grpbk = -1;
231	}
232	}
233
234	void makecompositiontable_p( short table, int pointt )
235	{
236	int point;
237
238	while( ( point = *pointt++ ) != END_OF_VEC )
239	table[point]++;
240	}
241
242
243	void makepointtable_nuc( int pointt, int n )
244	{
245	int point;
246	register int *p;
247
248	if( *n == -1 )
249	{
250	*pointt = -1;
251	return;
252	}
253
254	p = n;
255	point = n++ 1024;
256	point += n++ 256;
257	point += n++ 64;
258	point += n++ 16;
259	point += n++ 4;
260	point += *n++;
261	*pointt++ = point;
262
263	while( *n != END_OF_VEC )
264	{
265	point -= p++ 1024;
266	point *= 4;
267	point += *n++;
268	*pointt++ = point;
269	}
270	*pointt = END_OF_VEC;
271	}
272
273	void makepointtable( int pointt, int n )
274	{
275	int point;
276	register int *p;
277
278	if( *n == -1 )
279	{
280	*pointt = -1;
281	return;
282	}
283
284	p = n;
285	point = n++ 7776;
286	point += n++ 1296;
287	point += n++ 216;
288	point += n++ 36;
289	point += n++ 6;
290	point += *n++;
291	*pointt++ = point;
292
293	while( *n != END_OF_VEC )
294	{
295	point -= p++ 7776;
296	point *= 6;
297	point += *n++;
298	*pointt++ = point;
299	}
300	*pointt = END_OF_VEC;
301	}
302
303	static int commonsextet_p2( short table, int pointt )
304	{
305	int value = 0;
306	short tmp;
307	int point;
308	short *memo;
309	int *ct;
310	int *cp;
311
312	if( *pointt == -1 )
313	return( 0 );
314
315	memo = (short *)calloc( tsize, sizeof( short ) );
316	if( !memo ) ErrorExit( "Cannot allocate memo\n" );
317	ct = (int *)calloc( MIN( maxl, tsize )+1, sizeof( int ) ); // chuui!!
318	if( !ct ) ErrorExit( "Cannot allocate memo\n" );
319
320	cp = ct;
321	while( ( point = *pointt++ ) != END_OF_VEC )
322	{
323	tmp = memo[point]++;
324	if( tmp < table[point] )
325	value++;
326	if( tmp == 0 ) *cp++ = point;
327	}
328	*cp = END_OF_VEC;
329
330	cp = ct;
331	while( *cp != END_OF_VEC )
332	memo[*cp++] = 0;
333
334	free( memo );
335	free( ct );
336	return( value );
337	}
338
339	static int lastalign( char seq1, char seq2 )
340	{
341	FILE *lfp;
342	char command[10000];
343	int i, j, k;
344	int klim;
345	char kd[1000];
346	int msize = 10; // tekitou
347
348	lfp = fopen( "_scoringmatrixforlast", "w" ); // mochiron muda
349	if( !lfp )
350	{
351	fprintf( stderr, "Cannot open _scoringmatrixforlast" );
352	exit( 1 );
353	}
354	fprintf( lfp, " " );
355	for( j=0; j<4; j++ ) fprintf( lfp, " %c ", amino[j] );
356	fprintf( lfp, "\n" );
357	for( i=0; i<4; i++ )
358	{
359	fprintf( lfp, "%c ", amino[i] );
360	for( j=0; j<4; j++ ) fprintf( lfp, " %d ", n_dis[i][j] );
361	fprintf( lfp, "\n" );
362	}
363	fclose( lfp );
364
365	sprintf( command, "_dbd" );
366	lfp = fopen( command, "w" );
367	if( !lfp )
368	{
369	fprintf( stderr, "Cannot open _dbd" );
370	exit( 1 );
371	}
372	fprintf( lfp, ">0\n%s\n", seq1 );
373
374	fclose( lfp );
375	sprintf( command, "%s/lastdb _dbd _dbd", whereispairalign );
376	system( command );
377
378	sprintf( command, "_q", k );
379	lfp = fopen( command, "w" );
380	if( !lfp )
381	{
382	fprintf( stderr, "Cannot open %s", command );
383	exit( 1 );
384	}
385	fprintf( lfp, ">0\n%s\n", seq2 );
386	fclose( lfp );
387
388	sprintf( command, "%s/lastal -m %d -e %d -f 0 -s 1 -p _scoringmatrixforlast -a %d -b %d _db%sd _q > _lastres", whereispairalign, msize, laste, -penalty, -penalty_ex, kd );
389	if( system( command ) ) exit( 1 );
390
391
392	exit( 1 );
393	}
394
395	static void directionthread( void arg )
396	{
397	thread_arg_t targ = (thread_arg_t )arg;
398	int iend = targ->iend;
399	char **seq = targ->seq;
400	char *tmpseq = targ->tmpseq;
401	int *res = targ->res;
402	int **spointt = targ->spointt;
403	short *table1 = targ->table1;
404	int iq = targ->iq;
405	#ifdef enablemultithread
406	int thread_no = targ->thread_no;
407	int *jshare = targ->jshare;
408	#endif
409	int j;
410	char mseq1, mseq2;
411
412
413	if( dodp ) // nakuserukamo
414	{
415	mseq1 = AllocateCharMtx( 1, 0 );
416	mseq2 = AllocateCharMtx( 1, 0 );
417	}
418
419	j = -1;
420	while( 1 )
421	{
422	#ifdef enablemultithread
423	if( nthread )
424	{
425	pthread_mutex_lock( targ->mutex_counter );
426	j = *jshare;
427	if( j == iend )
428	{
429	fprintf( stderr, "\r %d / %d (thread %d) \r", iq, njob, thread_no );
430	pthread_mutex_unlock( targ->mutex_counter );
431	break;
432	}
433	++(*jshare);
434	pthread_mutex_unlock( targ->mutex_counter );
435	}
436	else
437	#endif
438	{
439	j++;
440	if( j == iend )
441	{
442	fprintf( stderr, "\r %d / %d \r", iq, njob );
443	break;
444	}
445	}
446
447
448	if( dodp == 1 )
449	{
450	// strcpy( mseq1[0], tmpseq );
451	// strcpy( mseq2[0], seq[j] );
452	mseq1[0] = tmpseq;
453	mseq2[0] = seq[j];
454	// res[j] = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, 0 );
455	res[j] = (int)L__align11_noalign( mseq1, mseq2 );
456	}
457	else if( dodp == 2 )
458	{
459	res[j] = lastalign( tmpseq, seq[j] );
460	}
461	else
462	res[j] = commonsextet_p2( table1, spointt[j] );
463	}
464	if( dodp ) // nakuserukamo
465	{
466	free( mseq1 );
467	free( mseq2 );
468	// G__align11_noalign( NULL, 0, 0, NULL, NULL, 0 );
469	L__align11_noalign( NULL, NULL );
470	}
471	// else
472	// if( nthread ) // inthread == 0 no toki free suru to, error. nazeda
473	// commonsextet_p( NULL, NULL );
474	return( NULL );
475	}
476
477	int main( int argc, char *argv[] )
478	{
479	static int *nlen;
480	static int *nogaplen;
481	static char name, seq;
482	int i, j, istart, iend;
483	FILE *infp;
484	// FILE *adfp;
485	char c;
486
487	int *grpseq;
488	char tmpseq, revseq;
489	int pointt, pointt_rev, **spointt;
490	float res_forward, res_reverse, res_max;
491	int ires, mres, mres2;
492	int *res;
493	static short table1, table1_rev;
494	static char mseq1f, mseq1r, **mseq2;
495
496	arguments( argc, argv );
497	#ifndef enablemultithread
498	nthread = 0;
499	#endif
500
501	if( inputfile )
502	{
503	infp = fopen( inputfile, "r" );
504	if( !infp )
505	{
506	fprintf( stderr, "Cannot open %s\n", inputfile );
507	exit( 1 );
508	}
509	}
510	else
511	infp = stdin;
512
513	getnumlen( infp );
514	rewind( infp );
515
516	if( alg == 'a' )
517	{
518	if( nlenmax < 10000 )
519	alg = 'G';
520	else
521	alg = 'S';
522	}
523
524	seq = AllocateCharMtx( njob, nlenmax*1+1 );
525
526	#if 0
527	Read( name, nlen, seq );
528	readData( infp, name, nlen, seq );
529	#else
530	name = AllocateCharMtx( njob, B+1 );
531	nlen = AllocateIntVec( njob );
532	nogaplen = AllocateIntVec( njob );
533	readData_pointer( infp, name, nlen, seq );
534	fclose( infp );
535
536	if( dorp != 'd' )
537	{
538	fprintf( stderr, "Not necessary!\n" );
539	for( i=0; i<njob; i++ )
540	fprintf( stdout, "_F_%-10.10s\n", name[i]+1 );
541	exit( 1 );
542	}
543	#endif
544
545	constants( njob, seq );
546
547
548	#if 0
549	fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
550	#endif
551
552	initSignalSM();
553
554	initFiles();
555
556	c = seqcheck( seq );
557	if( c )
558	{
559	fprintf( stderr, "Illegal character %c\n", c );
560	exit( 1 );
561	}
562
563	fprintf( stderr, "\n" );
564	if( alg == 'G' ) // do to the first sequence
565	{
566	mseq1f = AllocateCharMtx( 1, nlenmax+nlenmax );
567	mseq1r = AllocateCharMtx( 1, nlenmax+nlenmax );
568	mseq2 = AllocateCharMtx( 1, nlenmax+nlenmax );
569	tmpseq = AllocateCharVec( MAX( nlenmax, B ) +1 );
570
571	gappick0( mseq1f[0], seq[0] );
572	sreverse( mseq1r[0], mseq1f[0] );
573	strcpy( seq[0], mseq1f[0] );
574
575	if( nadd )
576	istart = njob - nadd;
577	else
578	istart = 1;
579
580	fprintf( stderr, "\n" );
581
582	for( i=0; i<istart; i++ )
583	{
584	gappick0( tmpseq, seq[i] );
585	strcpy( seq[i], tmpseq );
586	strcpy( tmpseq, name[i] );
587	strcpy( name[i], "_F_" );
588	strncpy( name[i]+3, tmpseq+1, 10 );
589	name[i][13] = 0;
590	}
591	for( i=istart; i<njob; i++ )
592	{
593	gappick0( mseq2[0], seq[i] );
594
595	// res_forward = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1f, mseq2, 0 );
596	// res_reverse = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1r, mseq2, 0 );
597
598	res_forward = L__align11_noalign( mseq1f, mseq2 );
599	res_reverse = L__align11_noalign( mseq1r, mseq2 );
600	#if 0
601
602	strcpy( mseq2[0], seq[i] );
603	strcpy( mseq1f[0], seq[0] );
604	res_forward = G__align11( mseq1f, mseq2, nlenmax*2, 0, 0 );
605	fprintf( stdout, "%s\n", mseq1f[0] );
606	fprintf( stdout, "%s\n", mseq2[0] );
607
608	strcpy( mseq2[0], seq[i] );
609	sreverse( mseq1r[0], seq[0] );
610	res_reverse = G__align11( mseq1r, mseq2, nlenmax*2, 0, 0 );
611	fprintf( stdout, "%s\n", mseq1r[0] );
612	fprintf( stdout, "%s\n", mseq2[0] );
613	#endif
614
615	// fprintf( stdout, "\nscore_for(%d,%d) = %f\n", 0, i, res_forward );
616	// fprintf( stdout, "score_rev(%d,%d) = %f\n", 0, i, res_reverse );
617	res_max = MAX(res_reverse,res_forward);
618	if( (res_reverse-res_forward)/res_max > thresholdtorev ) // tekitou
619	{
620	// fprintf( stderr, "REVERSE!!!\n" );
621	sreverse( seq[i], mseq2[0] );
622
623	strcpy( tmpseq, name[i] );
624	strcpy( name[i], "_R_" );
625	strncpy( name[i]+3, tmpseq+1, 10 );
626	name[i][13] = 0;
627	}
628	else
629	{
630	strcpy( seq[i], mseq2[0] );
631
632	strcpy( tmpseq, name[i] );
633	strcpy( name[i], "_F_" );
634	strncpy( name[i]+3, tmpseq+1, 10 );
635	name[i][13] = 0;
636	}
637	}
638	FreeCharMtx( mseq1f );
639	FreeCharMtx( mseq1r );
640	FreeCharMtx( mseq2 );
641	free( tmpseq );
642	}
643	else if( alg == 'm' )
644	{
645	if( dodp ) // nakuserukamo
646	{
647	mseq1f = AllocateCharMtx( 1, nlenmax+1);
648	mseq1r = AllocateCharMtx( 1, nlenmax+1 );
649	mseq2 = AllocateCharMtx( 1, nlenmax+1 );
650	}
651	else
652	{
653	spointt = AllocateIntMtx( njob, 0 );
654	pointt = AllocateIntMtx( njob, nlenmax+1 );
655	pointt_rev = AllocateIntMtx( njob, nlenmax+1 );
656	}
657	tmpseq = AllocateCharVec( MAX( nlenmax, B ) +1 );
658	revseq = AllocateCharVec( nlenmax+1 );
659	grpseq = AllocateIntVec( nlenmax+1 );
660	res = AllocateIntVec( njob );
661	if( dorp == 'd' ) tsize = (int)pow( 4, 6 );
662	else tsize = (int)pow( 6, 6 ); // iranai
663
664	maxl = 0;
665	for( i=0; i<njob; i++ )
666	{
667	gappick0( tmpseq, seq[i] );
668	nogaplen[i] = strlen( tmpseq );
669	if( nogaplen[i] > maxl ) maxl = nogaplen[i];
670	}
671
672	if( nadd )
673	iend = njob - nadd;
674	else
675	iend = 1;
676
677	for( i=0; i<iend; i++ )
678	{
679	// fprintf( stdout, "%d, SKIP\n", i );
680	gappick0( tmpseq, seq[i] );
681	strcpy( seq[i], tmpseq );
682	// if( !nadd ) strcpy( seq[i], tmpseq ); // seq ha tsukawanaikara ii.
683
684	if( !dodp )
685	{
686	seq_grp_nuc( grpseq, tmpseq );
687	makepointtable_nuc( pointt[i], grpseq );
688	spointt[i] = pointt[i];
689	}
690
691	strcpy( tmpseq, name[i] );
692	strcpy( name[i], "_F_" );
693	strncpy( name[i]+3, tmpseq+1, 10 );
694	name[i][13] = 0;
695	}
696
697	if( nadd )
698	istart = njob - nadd;
699	else
700	istart = 1;
701
702	fprintf( stderr, "\n" );
703	for( i=istart; i<njob; i++ )
704	{
705	// fprintf( stderr, "\r %d / %d ", i, njob );
706	gappick0( tmpseq, seq[i] );
707	strcpy( seq[i], tmpseq );
708	sreverse( revseq, tmpseq );
709
710	if( !dodp )
711	{
712	table1 = (short *)calloc( tsize, sizeof( short ) );
713	if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
714	table1_rev = (short *)calloc( tsize, sizeof( short ) );
715	if( !table1_rev ) ErrorExit( "Cannot allocate table1_rev\n" );
716	seq_grp_nuc( grpseq, tmpseq );
717	makepointtable_nuc( pointt[i], grpseq );
718	makecompositiontable_p( table1, pointt[i] );
719	seq_grp_nuc( grpseq, revseq );
720	makepointtable_nuc( pointt_rev[i], grpseq );
721	makecompositiontable_p( table1_rev, pointt_rev[i] );
722	}
723
724	if( nadd && addfragment )
725	iend = njob-nadd;
726	else
727	iend = i;
728
729
730	#ifdef enablemultithread
731	if( nthread )
732	{
733	pthread_t *handle;
734	pthread_mutex_t mutex_counter;
735	thread_arg_t *targ;
736	int *jsharept;
737
738	targ = calloc( nthread, sizeof( thread_arg_t ) );
739	handle = calloc( nthread, sizeof( pthread_t ) );
740	pthread_mutex_init( &mutex_counter, NULL );
741	jsharept = calloc( 1, sizeof(int) );
742	*jsharept = 0;
743
744	for( j=0; j<nthread; j++ )
745	{
746	targ[j].iend = iend;
747	targ[j].seq = seq;
748	targ[j].tmpseq = tmpseq;
749	targ[j].res = res;
750	targ[j].spointt = spointt;
751	targ[j].table1 = table1;
752	targ[j].jshare = jsharept;
753	targ[j].iq = i;
754	targ[j].mutex_counter = &mutex_counter;
755	targ[j].thread_no = j;
756	pthread_create( handle+j, NULL, directionthread, (void *)(targ+j) );
757	}
758	for( j=0; j<nthread; j++ ) pthread_join( handle[j], NULL );
759	pthread_mutex_destroy( &mutex_counter );
760	free( handle );
761	free( targ );
762	free( jsharept );
763	}
764	else
765	#endif
766	{
767	thread_arg_t *targ;
768	targ = calloc( 1, sizeof( thread_arg_t ) );
769	targ[0].iend = iend;
770	targ[0].seq = seq;
771	targ[0].tmpseq = tmpseq;
772	targ[0].res = res;
773	targ[0].spointt = spointt;
774	targ[0].table1 = table1;
775	targ[0].iq = i;
776	directionthread( targ );
777	free( targ );
778	}
779
780
781	mres = mres2 = 0;
782	for( j=0; j<iend; j++ )
783	{
784	ires = res[j];
785	// fprintf( stdout, "ires (%d,%d) = %d\n", i, j, ires );
786	// fflush( stdout );
787	if( ires>mres2 )
788	{
789	if( ires>mres )
790	{
791	mres2 = mres;
792	mres = ires;
793	}
794	else
795	mres2 = ires;
796	}
797	}
798	res_forward = (float)( mres + mres2 ) / 2;
799
800	#ifdef enablemultithread
801	if( nthread )
802	{
803	pthread_t *handle;
804	pthread_mutex_t mutex_counter;
805	thread_arg_t *targ;
806	int *jsharept;
807
808	targ = calloc( nthread, sizeof( thread_arg_t ) );
809	handle = calloc( nthread, sizeof( pthread_t ) );
810	pthread_mutex_init( &mutex_counter, NULL );
811	jsharept = calloc( 1, sizeof(int) );
812	*jsharept = 0;
813
814	for( j=0; j<nthread; j++ )
815	{
816	targ[j].iend = iend;
817	targ[j].seq = seq;
818	targ[j].tmpseq = revseq;
819	targ[j].res = res;
820	targ[j].spointt = spointt;
821	targ[j].table1 = table1_rev;
822	targ[j].jshare = jsharept;
823	targ[j].iq = i;
824	targ[j].mutex_counter = &mutex_counter;
825	targ[j].thread_no = j;
826	pthread_create( handle+j, NULL, directionthread, (void *)(targ+j) );
827	}
828	for( j=0; j<nthread; j++ ) pthread_join( handle[j], NULL );
829	pthread_mutex_destroy( &mutex_counter );
830	free( handle );
831	free( targ );
832	free( jsharept );
833	}
834	else
835	#endif
836	{
837	thread_arg_t *targ;
838	targ = calloc( 1, sizeof( thread_arg_t ) );
839	targ[0].iend = iend;
840	targ[0].seq = seq;
841	targ[0].tmpseq = revseq;
842	targ[0].res = res;
843	targ[0].spointt = spointt;
844	targ[0].table1 = table1_rev;
845	targ[0].iq = i;
846	directionthread( targ );
847	free( targ );
848	}
849
850	mres = mres2 = 0;
851	for( j=0; j<iend; j++ )
852	{
853	ires = res[j];
854	if( ires>mres2 )
855	{
856	if( ires>mres )
857	{
858	mres2 = mres;
859	mres = ires;
860	}
861	else
862	mres2 = ires;
863	}
864	}
865	res_reverse = (float)( mres + mres2 ) / 2;
866
867	// fprintf( stdout, "\n" );
868	// fprintf( stdout, "score_for(%d,%d) = %f\n", 0, i, res_forward );
869	// fprintf( stdout, "score_rev(%d,%d) = %f\n", 0, i, res_reverse );
870	// fflush( stdout );
871	res_max = MAX(res_reverse,res_forward);
872	if( (res_reverse-res_forward)/res_max > thresholdtorev ) // tekitou
873
874	{
875	strcpy( seq[i], revseq );
876
877	strcpy( tmpseq, name[i] );
878	strcpy( name[i], "_R_" );
879	strncpy( name[i]+3, tmpseq+1, 10 );
880	name[i][13] = 0;
881	if( !dodp ) spointt[i] = pointt_rev[i];
882	}
883	else
884	{
885	strcpy( tmpseq, name[i] );
886	strcpy( name[i], "_F_" );
887	strncpy( name[i]+3, tmpseq+1, 10 );
888	name[i][13] = 0;
889	if( !dodp ) spointt[i] = pointt[i];
890	}
891
892	if( !dodp )
893	{
894	free( table1 );
895	free( table1_rev );
896	}
897	}
898
899	free( grpseq );
900	free( tmpseq );
901	free( revseq );
902	free( res );
903	if( dodp )
904	{
905	FreeCharMtx( mseq1f );
906	FreeCharMtx( mseq1r );
907	FreeCharMtx( mseq2 );
908	}
909	else
910	{
911	FreeIntMtx( pointt );
912	FreeIntMtx( pointt_rev );
913	free( spointt );
914	}
915	}
916	else
917	{
918	fprintf( stderr, "Unknown alg %c\n", alg );
919	exit( 1 );
920	}
921	// writeData_pointer( stdout, njob, name, nlen, seq );
922	for( i=0; i<njob; i++ )
923	{
924	// fprintf( stdout, ">%s\n", name[i] );
925	// fprintf( stdout, "%s\n", seq[i] );
926	fprintf( stdout, "%s\n", name[i] );
927	}
928
929	fprintf( stderr, "\n" );
930	SHOWVERSION;
931	return( 0 );
932	}
933

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

Download in other formats:

Original Format