1 /* Tree-dependent-iteration */
2 /* Devide to segments */
9 void arguments( int argc, char *argv[] )
24 divWinSize = 20; /* 70 */
39 alg = 'A'; /* chuui */
46 ppenalty = NOTSPECIFIED;
47 ppenalty_ex = NOTSPECIFIED;
48 poffset = NOTSPECIFIED;
49 kimuraR = NOTSPECIFIED;
52 fftWinSize = NOTSPECIFIED;
53 fftThreshold = NOTSPECIFIED;
57 while( --argc > 0 && (*++argv)[0] == '-' )
59 while ( c = *++argv[0] )
65 fprintf( stderr, "inputfile = %s\n", inputfile );
69 niter = atoi( *++argv );
70 fprintf( stderr, "niter = %d\n", niter );
74 ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
75 fprintf( stderr, "ppenalty = %d\n", ppenalty );
79 ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
80 fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
84 poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
85 fprintf( stderr, "poffset = %d\n", poffset );
89 kimuraR = atoi( *++argv );
90 fprintf( stderr, "kimuraR = %d\n", kimuraR );
94 nblosum = atoi( *++argv );
96 fprintf( stderr, "blosum %d\n", nblosum );
100 pamN = atoi( *++argv );
103 fprintf( stderr, "jtt %d\n", pamN );
107 pamN = atoi( *++argv );
110 fprintf( stderr, "tm %d\n", pamN );
114 fastathreshold = atof( *++argv );
116 fprintf( stderr, "fastathreshold %f\n", fastathreshold );
200 fftThreshold = atoi( *++argv );
204 fftWinSize = atoi( *++argv );
208 fprintf( stderr, "illegal option %c\n", c );
218 cut = atof( (*argv) );
223 fprintf( stderr, "options : Check source file!\n" );
227 if( alg == 'A' && weight == 0 )
228 ErrorExit( "ERROR : Algorithm A+ and un-weighted\n" );
233 int main( int argc, char *argv[] )
237 static char name[M][B], **seq, **aseq, **bseq;
238 static Segment *segment = NULL;
239 static int anchors[MAXSEG];
252 LocalHom **localhomtable;
254 arguments( argc, argv );
258 infp = fopen( inputfile, "r" );
261 fprintf( stderr, "Cannot open %s\n", inputfile );
269 PreRead( stdin, &njob, &nlenmax );
277 fprintf( stderr, "At least 2 sequences should be input!\n"
278 "Only %d sequence found.\n", njob );
284 segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
285 topol = AllocateIntCub( njob, 2, njob );
286 len = AllocateDoubleMtx( njob, 2 );
287 eff = AllocateDoubleMtx( njob, njob );
288 seq = AllocateCharMtx( njob, nlenmax*5+1 );
289 seq_g = AllocateCharMtx( njob, nlenmax*5+1 );
290 res_g = AllocateCharMtx( njob, nlenmax*5+1 );
291 aseq = AllocateCharMtx( njob, nlenmax*5+1 );
292 bseq = AllocateCharMtx( njob, nlenmax*5+1 );
293 alloclen = nlenmax * 5;
294 seq_g_bk = AllocateCharMtx( njob, 0 );
295 for( i=0; i<njob; i++ ) seq_g_bk[i] = seq_g[i];
299 localhomtable = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
300 for( i=0; i<njob; i++)
302 localhomtable[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
303 for( j=0; j<njob; j++)
305 localhomtable[i][j].start1 = -1;
306 localhomtable[i][j].end1 = -1;
307 localhomtable[i][j].start2 = -1;
308 localhomtable[i][j].end2 = -1;
309 localhomtable[i][j].overlapaa = -1.0;
310 localhomtable[i][j].opt = -1.0;
311 localhomtable[i][j].importance = -1.0;
312 localhomtable[i][j].next = NULL;
315 fprintf( stderr, "Loading 'hat3' ... " );
316 prep = fopen( "hat3", "r" );
317 if( prep == NULL ) ErrorExit( "Make hat3." );
318 readlocalhomtable( prep, njob, localhomtable );
320 // for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ )
321 // fprintf( stdout, "%d %d %d %d %d %d %d\n", i, j, localhomtable[i][j].opt, localhomtable[i][j].start1, localhomtable[i][j].end1, localhomtable[i][j].start2, localhomtable[i][j].end2 );
322 fprintf( stderr, "done.\n" );
324 fprintf( stderr, "Extending localhom ... " );
325 extendlocalhom( njob, localhomtable );
326 fprintf( stderr, "done.\n" );
331 Read( name, nlen, seq_g );
333 readData( infp, name, nlen, seq_g );
338 for( i=0; i<njob; i++ )
344 for( i=0; i<njob; i++ )
346 identity *= ( nlen[i] == nlen[0] );
350 fprintf( stderr, "Input pre-aligned data\n" );
353 constants( njob, seq_g );
356 fprintf( stderr, "penalty = %d\n", penalty );
357 fprintf( stderr, "penalty_ex = %d\n", penalty_ex );
358 fprintf( stderr, "offset = %d\n", offset );
368 writePre( njob, name, nlen, seq_g, 1 );
374 c = seqcheck( seq_g );
377 fprintf( stderr, "Illeagal character %c\n", c );
380 commongappick( njob, seq_g );
386 prep = fopen( "hat2", "r" );
387 if( !prep ) ErrorExit( "Make hat2." );
388 readhat2( prep, njob, name, eff );
391 for( i=0; i<njob-1; i++ )
393 for( j=i+1; j<njob; j++ )
395 printf( " %f", eff[i][j] );
400 if ( treemethod == 'x' )
401 veryfastsupg( njob, eff, topol, len );
402 else if( treemethod == 'n' )
403 nj( njob, eff, topol, len );
404 else if( treemethod == 's' )
405 spg( njob, eff, topol, len );
406 else if( treemethod == 'p' )
407 upg2( njob, eff, topol, len );
408 else ErrorExit( "Incorrect treemethod.\n" );
411 printf( "utree = %d\n", utree );
414 if( ( !utree && kobetsubunkatsu ) || constraint || !bunkatsu )
418 anchors[1] =strlen( seq_g[0] );
423 nseg = searchAnchors( njob, seq_g, segment );
425 fprintf( stderr, "nseg = %d\n", nseg );
426 fprintf( stderr, "seq_g[0] = %s\n", seq_g[0] );
427 fprintf( stderr, "nlenmax = %d\n", nlenmax );
428 fprintf( stderr, "len = %d\n", strlen( seq_g[0] ) );
432 for( i=0; i<nseg; i++ ) anchors[i+1] = segment[i].center;
433 anchors[nseg+1] = strlen( seq_g[0] );
437 for( i=0; i<nseg; i++ )
438 fprintf( stderr, "anchor[%d] = %d\n", i, anchors[i] );
442 for( i=0; i<njob; i++ ) res_g[i][0] = 0;
444 for( iseg=0; iseg<nseg-1; iseg++ )
446 int tmplen = anchors[iseg+1]-anchors[iseg];
447 int pos = strlen( res_g[0] );
448 for( j=0; j<njob; j++ )
450 strncpy( seq[j], seq_g[j], tmplen );
455 fprintf( stderr, "Segment %3d/%3d %4d-%4d\n", iseg+1, nseg-1, pos+1, pos+1+tmplen );
456 fprintf( trap_g, "Segment %3d/%3d %4d-%4d\n", iseg+1, nseg-1, pos+1, pos+1+tmplen );
459 returnvalue = TreeDependentIteration( njob, name, nlen, seq, bseq, topol, len, alloclen, localhomtable );
461 for( i=0; i<njob; i++ )
462 strcat( res_g[i], bseq[i] );
464 FreeCharMtx( seq_g_bk );
466 FreeDoubleMtx( len );
467 FreeDoubleMtx( eff );
468 fprintf( stderr, "constraint = %d\n", constraint );
469 if( constraint ) FreeLocalHomTable( localhomtable, njob );
472 Write( stdout, njob, name, nlen, bseq );
475 fprintf( stderr, "done\n" );
476 fprintf( trap_g, "done\n" );
481 writePre( njob, name, nlen, res_g, 1 );
483 writeData( stdout, njob, name, nlen, res_g, 1 );
492 signed int main( int argc, char *argv[] )
499 gets( b ); njob = atoi( b );
503 if( strstr( b, "ayhoff" ) ) scoremtx = 1;
504 else if( strstr( b, "dna" ) || strstr( b, "DNA" ) ) scoremtx = -1;
505 else if( strstr( b, "M-Y" ) || strstr( b, "iyata" ) ) scoremtx = 2;
508 if( strstr( b, "constraint" ) ) cnst = 1;
515 if( !strncmp( b, a, 1 ) )
517 gets( b ); nlen[i] = atoi( b );
518 if( nlen[i] > nlenmax ) nlenmax = nlen[i];
522 if( nlenmax > N || njob > M )
524 fprintf( stderr, "ERROR in main\n" );
531 value = main1( nlen, argc, argv );