Context Navigation

rnatest.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: 9.8 KB

Line
1	#include "mltaln.h"
2
3	#define DEBUG 0
4	#define IODEBUG 0
5	#define SCOREOUT 0
6
7
8
9	void arguments( int argc, char *argv[] )
10	{
11	int c;
12
13	inputfile = NULL;
14	fftkeika = 0;
15	pslocal = -1000.0;
16	constraint = 0;
17	nblosum = 62;
18	fmodel = 0;
19	calledByXced = 0;
20	devide = 0;
21	use_fft = 0;
22	fftscore = 1;
23	fftRepeatStop = 0;
24	fftNoAnchStop = 0;
25	weight = 3;
26	utree = 1;
27	tbutree = 1;
28	refine = 0;
29	check = 1;
30	cut = 0.0;
31	disp = 0;
32	outgap = 1;
33	alg = 'A';
34	mix = 0;
35	tbitr = 0;
36	scmtd = 5;
37	tbweight = 0;
38	tbrweight = 3;
39	checkC = 0;
40	treemethod = 'x';
41	contin = 0;
42	scoremtx = 1;
43	kobetsubunkatsu = 0;
44	divpairscore = 0;
45	dorp = NOTSPECIFIED;
46	ppenalty = NOTSPECIFIED;
47	ppenalty_OP = NOTSPECIFIED;
48	ppenalty_ex = NOTSPECIFIED;
49	ppenalty_EX = NOTSPECIFIED;
50	poffset = NOTSPECIFIED;
51	kimuraR = NOTSPECIFIED;
52	pamN = NOTSPECIFIED;
53	geta2 = GETA2;
54	fftWinSize = NOTSPECIFIED;
55	fftThreshold = NOTSPECIFIED;
56	RNAppenalty = NOTSPECIFIED;
57	RNApthr = NOTSPECIFIED;
58
59	while( --argc > 0 && (*++argv)[0] == '-' )
60	{
61	while ( ( c = *++argv[0] ) )
62	{
63	switch( c )
64	{
65	case 'i':
66	inputfile = *++argv;
67	fprintf( stderr, "inputfile = %s\n", inputfile );
68	--argc;
69	goto nextoption;
70	case 'o':
71	RNAppenalty = (int)( atof( ++argv ) 1000 - 0.5 );
72	--argc;
73	goto nextoption;
74	case 'f':
75	ppenalty = (int)( atof( ++argv ) 1000 - 0.5 );
76	--argc;
77	goto nextoption;
78	case 'g':
79	ppenalty_ex = (int)( atof( ++argv ) 1000 - 0.5 );
80	--argc;
81	goto nextoption;
82	case 'O':
83	ppenalty_OP = (int)( atof( ++argv ) 1000 - 0.5 );
84	--argc;
85	goto nextoption;
86	case 'E':
87	ppenalty_EX = (int)( atof( ++argv ) 1000 - 0.5 );
88	--argc;
89	goto nextoption;
90	case 'h':
91	poffset = (int)( atof( ++argv ) 1000 - 0.5 );
92	--argc;
93	goto nextoption;
94	case 'k':
95	kimuraR = atoi( *++argv );
96	// fprintf( stderr, "kimuraR = %d\n", kimuraR );
97	--argc;
98	goto nextoption;
99	case 'b':
100	nblosum = atoi( *++argv );
101	scoremtx = 1;
102	fprintf( stderr, "blosum %d\n", nblosum );
103	--argc;
104	goto nextoption;
105	case 'j':
106	pamN = atoi( *++argv );
107	scoremtx = 0;
108	TMorJTT = JTT;
109	fprintf( stderr, "jtt %d\n", pamN );
110	--argc;
111	goto nextoption;
112	case 'm':
113	pamN = atoi( *++argv );
114	scoremtx = 0;
115	TMorJTT = TM;
116	fprintf( stderr, "TM %d\n", pamN );
117	--argc;
118	goto nextoption;
119	case 'l':
120	ppslocal = (int)( atof( ++argv ) 1000 + 0.5 );
121	pslocal = (int)( 600.0 / 1000.0 * ppslocal + 0.5);
122	// fprintf( stderr, "ppslocal = %d\n", ppslocal );
123	// fprintf( stderr, "pslocal = %d\n", pslocal );
124	--argc;
125	goto nextoption;
126	#if 1
127	case 'a':
128	fmodel = 1;
129	break;
130	#endif
131	case 'r':
132	fmodel = -1;
133	break;
134	case 'D':
135	dorp = 'd';
136	break;
137	case 'P':
138	dorp = 'p';
139	break;
140	#if 0
141	case 'e':
142	fftscore = 0;
143	break;
144	case 'O':
145	fftNoAnchStop = 1;
146	break;
147	case 'R':
148	fftRepeatStop = 1;
149	break;
150	#endif
151	case 'Q':
152	calledByXced = 1;
153	break;
154	case 's':
155	treemethod = 's';
156	break;
157	case 'x':
158	disp = 1;
159	break;
160	case 'p':
161	treemethod = 'p';
162	break;
163	#if 0
164	case 'a':
165	alg = 'a';
166	break;
167	#endif
168	case 'S':
169	alg = 'S';
170	break;
171	case 'L':
172	alg = 'L';
173	break;
174	case 'M':
175	alg = 'M';
176	break;
177	case 'R':
178	alg = 'R';
179	break;
180	case 'N':
181	alg = 'N';
182	break;
183	case 'A':
184	alg = 'A';
185	break;
186	case 'V':
187	alg = 'V';
188	break;
189	case 'C':
190	alg = 'C';
191	break;
192	case 'F':
193	use_fft = 1;
194	break;
195	case 'v':
196	tbrweight = 3;
197	break;
198	case 'd':
199	divpairscore = 1;
200	break;
201	/* Modified 01/08/27, default: user tree */
202	case 'J':
203	tbutree = 0;
204	break;
205	/* modification end. */
206	case 'z':
207	fftThreshold = atoi( *++argv );
208	--argc;
209	goto nextoption;
210	case 'w':
211	fftWinSize = atoi( *++argv );
212	--argc;
213	goto nextoption;
214	case 'Z':
215	checkC = 1;
216	break;
217	default:
218	fprintf( stderr, "illegal option %c\n", c );
219	argc = 0;
220	break;
221	}
222	}
223	nextoption:
224	;
225	}
226	if( argc == 1 )
227	{
228	cut = atof( (*argv) );
229	argc--;
230	}
231	if( argc != 0 )
232	{
233	fprintf( stderr, "options: Check source file !\n" );
234	exit( 1 );
235	}
236	if( tbitr == 1 && outgap == 0 )
237	{
238	fprintf( stderr, "conflicting options : o, m or u\n" );
239	exit( 1 );
240	}
241	if( alg == 'C' && outgap == 0 )
242	{
243	fprintf( stderr, "conflicting options : C, o\n" );
244	exit( 1 );
245	}
246	}
247
248	int countamino( char *s, int end )
249	{
250	int val = 0;
251	while( end-- )
252	if( *s++ != '-' ) val++;
253	return( val );
254	}
255
256
257	static void WriteOptions( FILE *fp )
258	{
259
260	if( dorp == 'd' ) fprintf( fp, "DNA\n" );
261	else
262	{
263	if ( scoremtx == 0 ) fprintf( fp, "JTT %dPAM\n", pamN );
264	else if( scoremtx == 1 ) fprintf( fp, "BLOSUM %d\n", nblosum );
265	else if( scoremtx == 2 ) fprintf( fp, "M-Y\n" );
266	}
267	fprintf( stderr, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
268	if( use_fft ) fprintf( fp, "FFT on\n" );
269
270	fprintf( fp, "tree-base method\n" );
271	if( tbrweight == 0 ) fprintf( fp, "unweighted\n" );
272	else if( tbrweight == 3 ) fprintf( fp, "clustalw-like weighting\n" );
273	if( tbitr \|\| tbweight )
274	{
275	fprintf( fp, "iterate at each step\n" );
276	if( tbitr && tbrweight == 0 ) fprintf( fp, " unweighted\n" );
277	if( tbitr && tbrweight == 3 ) fprintf( fp, " reversely weighted\n" );
278	if( tbweight ) fprintf( fp, " weighted\n" );
279	fprintf( fp, "\n" );
280	}
281
282	fprintf( fp, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
283
284	if( alg == 'a' )
285	fprintf( fp, "Algorithm A\n" );
286	else if( alg == 'A' )
287	fprintf( fp, "Algorithm A+\n" );
288	else if( alg == 'S' )
289	fprintf( fp, "Apgorithm S\n" );
290	else if( alg == 'C' )
291	fprintf( fp, "Apgorithm A+/C\n" );
292	else
293	fprintf( fp, "Unknown algorithm\n" );
294
295	if( treemethod == 'x' )
296	fprintf( fp, "Tree = UPGMA (3).\n" );
297	else if( treemethod == 's' )
298	fprintf( fp, "Tree = UPGMA (2).\n" );
299	else if( treemethod == 'p' )
300	fprintf( fp, "Tree = UPGMA (1).\n" );
301	else
302	fprintf( fp, "Unknown tree.\n" );
303
304	if( use_fft )
305	{
306	fprintf( fp, "FFT on\n" );
307	if( dorp == 'd' )
308	fprintf( fp, "Basis : 4 nucleotides\n" );
309	else
310	{
311	if( fftscore )
312	fprintf( fp, "Basis : Polarity and Volume\n" );
313	else
314	fprintf( fp, "Basis : 20 amino acids\n" );
315	}
316	fprintf( fp, "Threshold of anchors = %d%%\n", fftThreshold );
317	fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
318	}
319	else
320	fprintf( fp, "FFT off\n" );
321	fflush( fp );
322	}
323
324
325	int main( int argc, char *argv[] )
326	{
327	static int *nlen;
328	static char name, seq, **useq;
329	static char mseq1, mseq2;
330	static char **aseq;
331	static char **bseq;
332	static double *eff;
333	int i;
334	FILE *infp;
335	char c;
336	int alloclen;
337	RNApair **pair1;
338	RNApair **pair2;
339	float **map;
340
341	arguments( argc, argv );
342
343	if( inputfile )
344	{
345	infp = fopen( inputfile, "r" );
346	if( !infp )
347	{
348	fprintf( stderr, "Cannot open %s\n", inputfile );
349	exit( 1 );
350	}
351	}
352	else
353	infp = stdin;
354
355	getnumlen( infp );
356	rewind( infp );
357
358	if( njob > M )
359	{
360	fprintf( stderr, "The number of sequences must be < %d\n", M );
361	fprintf( stderr, "Please try the splittbfast program for such large data.\n" );
362	exit( 1 );
363	}
364
365	name = AllocateCharMtx( njob, B+1 );
366	nlen = AllocateIntVec( njob );
367
368	seq = AllocateCharMtx( njob, nlenmax*5+1 );
369	useq = AllocateCharMtx( njob, nlenmax*5+1 );
370	aseq = AllocateCharMtx( njob, nlenmax*5+1 );
371	bseq = AllocateCharMtx( njob, nlenmax*5+1 );
372	mseq1 = AllocateCharMtx( njob, 0 );
373	mseq2 = AllocateCharMtx( njob, 0 );
374	alloclen = nlenmax*5;
375
376	pair1 = calloc( nlenmax5+1, sizeof( RNApair ) );
377	pair2 = calloc( nlenmax5+1, sizeof( RNApair ) );
378	map = AllocateFloatMtx( nlenmax+1, nlenmax );
379
380	eff = AllocateDoubleVec( njob );
381
382	readData_pointer( infp, name, nlen, seq );
383	fclose( infp );
384
385	for( i=0; i<njob; i++ ) strcpy( useq[i], seq[i] );
386
387	constants( njob, seq );
388
389	#if 0
390	fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
391	#endif
392
393	initSignalSM();
394
395	initFiles();
396
397	WriteOptions( trap_g );
398
399	c = seqcheck( seq );
400	if( c )
401	{
402	fprintf( stderr, "Illeagal character %c\n", c );
403	exit( 1 );
404	}
405
406	// writePre( njob, name, nlen, seq, 0 );
407
408	for( i=0; i<njob; i++ ) eff[i] = 1.0;
409
410	// for( i=0; i<njob; i++ ) gappick0( bseq[i], seq[i] );
411
412
413	fprintf( stderr, "folding group1\n" );
414	// rnalocal( seq, useq, eff, eff, njob, njob, alloclen, pair1 );
415	rnaalifoldcall( seq, njob, pair1 );
416	exit( 1 );
417
418	#if 0
419	fprintf( stderr, "folding group1\n" );
420	rnalocal( seq+1, useq+1, eff+1, eff+1, 1, 1, alloclen, pair2 );
421	fprintf( stderr, "aligning 1 and 2, phase 1\n" );
422	Lalignmm_hmout( seq, seq+1, eff, eff+1, 1, 1, alloclen, NULL, NULL, NULL, NULL, map );
423
424
425	#if 0
426	lgth1 = strlen( seq[0] );
427	for( i=0; i<lgth1; i++ )
428	{
429	fprintf( stderr, "\n" );
430	if( pair1[i].pos == -1 ) continue;
431	fprintf( stderr, "i=%d (%c):%d", i, seq[0][i], pair1[i].pos );
432	if( map12[pair1[i].pos].pos == -1 ) continue;
433	fprintf( stderr, "%c -> %c ", seq[0][pair1[i].pos], seq[1][map12[pair1[i].pos].pos] );
434	if( pair2[map12[pair1[i].pos].pos].pos == -1 ) continue;
435	fprintf( stderr, "%d:%d (%c)", map12[pair1[i].pos].pos, pair2[map12[pair1[i].pos].pos].pos, seq[1][pair2[map12[pair1[i].pos].pos].pos] );
436	}
437	#endif
438
439
440	exit( 1 );
441
442
443	pairalign( name, nlen, bseq, aseq, mseq1, mseq2, eff, alloclen );
444	fprintf( trap_g, "done.\n" );
445	#if DEBUG
446	fprintf( stderr, "closing trap_g\n" );
447	#endif
448	fclose( trap_g );
449
450	// writePre( njob, name, nlen, aseq, !contin );
451	#if 0
452	writeData( stdout, njob, name, nlen, aseq );
453	#endif
454	#if IODEBUG
455	fprintf( stderr, "OSHIMAI\n" );
456	#endif
457	SHOWVERSION;
458	return( 0 );
459	#endif
460	}

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

Context Navigation

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

Download in other formats: