1 /* Tree-dependent-iteration */
2 /* Devide to segments */
12 void arguments( int argc, char *argv[] )
21 parallelizationstrategy = BAATARI1;
38 divWinSize = 20; /* 70 */
52 use_fft = 0; // CHUUI dochira demo mafft.tmpl deha F
54 alg = 'A'; /* chuui */
61 ppenalty = NOTSPECIFIED;
62 ppenalty_ex = NOTSPECIFIED;
63 poffset = NOTSPECIFIED;
64 kimuraR = NOTSPECIFIED;
67 fftWinSize = NOTSPECIFIED;
68 fftThreshold = NOTSPECIFIED;
69 RNAppenalty = NOTSPECIFIED;
70 RNAppenalty_ex = NOTSPECIFIED;
71 RNApthr = NOTSPECIFIED;
73 consweight_multi = 1.0;
76 while( --argc > 0 && (*++argv)[0] == '-' )
78 while ( (c = *++argv[0]) )
84 fprintf( stderr, "inputfile = %s\n", inputfile );
88 niter = atoi( *++argv );
89 fprintf( stderr, "niter = %d\n", niter );
93 RNApthr = (int)( atof( *++argv ) * 1000 - 0.5 );
97 RNAppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
101 ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
102 // fprintf( stderr, "ppenalty = %d\n", ppenalty );
106 ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
107 // fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
111 poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
112 fprintf( stderr, "poffset = %d\n", poffset );
116 kimuraR = atoi( *++argv );
117 fprintf( stderr, "kappa = %d\n", kimuraR );
121 nblosum = atoi( *++argv );
123 fprintf( stderr, "blosum %d / kimura 200\n", nblosum );
127 pamN = atoi( *++argv );
130 fprintf( stderr, "jtt/kimura %d\n", pamN );
134 pamN = atoi( *++argv );
137 fprintf( stderr, "tm %d\n", pamN );
141 fastathreshold = atof( *++argv );
146 consweight_rna = atof( *++argv );
151 consweight_multi = atof( *++argv );
155 nthread = atoi( *++argv );
156 fprintf( stderr, "nthread = %d\n", nthread );
160 randomseed = atoi( *++argv );
161 fprintf( stderr, "randomseed = %d\n", randomseed );
166 if( !strcmp( argkey, "BESTFIRST" ) ) parallelizationstrategy = BESTFIRST;
167 else if( !strcmp( argkey, "BAATARI0" ) ) parallelizationstrategy = BAATARI0;
168 else if( !strcmp( argkey, "BAATARI1" ) ) parallelizationstrategy = BAATARI1;
169 else if( !strcmp( argkey, "BAATARI2" ) ) parallelizationstrategy = BAATARI2;
172 fprintf( stderr, "Unknown parallelization strategy, %s\n", argkey );
291 fftThreshold = atoi( *++argv );
295 fftWinSize = atoi( *++argv );
299 fprintf( stderr, "illegal option %c\n", c );
309 cut = atof( (*argv) );
314 fprintf( stderr, "options : Check source file!\n" );
318 if( alg == 'A' && weight == 0 )
319 ErrorExit( "ERROR : Algorithm A+ and un-weighted\n" );
324 int main( int argc, char *argv[] )
328 static char **name, **seq, **aseq, **bseq;
329 static Segment *segment = NULL;
330 static int anchors[MAXSEG];
343 LocalHom **localhomtable = NULL; // by D.Mathog
344 RNApair ***singlerna;
346 static char **nogap1seq;
347 static char *kozoarivec;
350 arguments( argc, argv );
351 #ifndef enablemultithread
357 infp = fopen( inputfile, "r" );
360 fprintf( stderr, "Cannot open %s\n", inputfile );
368 PreRead( stdin, &njob, &nlenmax );
378 fprintf( stderr, "At least 2 sequences should be input!\n"
379 "Only %d sequence found.\n", njob );
385 segment = (Segment *)calloc( MAXSEG, sizeof( Segment ) );
386 // if( treemethod == 'X' || treemethod == 'E' || treemethod == 'q' )
387 topol = AllocateIntCub( njob, 2, njob );
388 len = AllocateDoubleMtx( njob, 2 );
389 eff = AllocateDoubleMtx( njob, njob );
390 seq = AllocateCharMtx( njob, nlenmax*9+1 );
391 name = AllocateCharMtx( njob, B+1 );
392 seq_g = AllocateCharMtx( njob, nlenmax*9+1 );
393 res_g = AllocateCharMtx( njob, nlenmax*9+1 );
394 aseq = AllocateCharMtx( njob, nlenmax*9+1 );
395 bseq = AllocateCharMtx( njob, nlenmax*9+1 );
396 nogap1seq = AllocateCharMtx( 1, nlenmax*1+1 );
397 alloclen = nlenmax * 9;
398 seq_g_bk = AllocateCharMtx( njob, 0 );
399 for( i=0; i<njob; i++ ) seq_g_bk[i] = seq_g[i];
400 kozoarivec = AllocateCharVec( njob );
404 localhomtable = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
405 for( i=0; i<njob; i++)
407 localhomtable[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
408 for( j=0; j<njob; j++)
410 localhomtable[i][j].start1 = -1;
411 localhomtable[i][j].end1 = -1;
412 localhomtable[i][j].start2 = -1;
413 localhomtable[i][j].end2 = -1;
414 localhomtable[i][j].overlapaa = -1.0;
415 localhomtable[i][j].opt = -1.0;
416 localhomtable[i][j].importance = -1.0;
417 localhomtable[i][j].next = NULL;
418 localhomtable[i][j].nokori = 0;
419 localhomtable[i][j].extended = -1;
420 localhomtable[i][j].last = localhomtable[i]+j;
421 localhomtable[i][j].korh = 'h';
424 fprintf( stderr, "Loading 'hat3' ... " );
426 prep = fopen( "hat3", "r" );
427 if( prep == NULL ) ErrorExit( "Make hat3." );
428 readlocalhomtable2( prep, njob, localhomtable, kozoarivec );
430 // for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ )
431 // 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 );
432 fprintf( stderr, "done.\n" );
435 fprintf( stderr, "Extending localhom ... " );
436 extendlocalhom( njob, localhomtable );
437 fprintf( stderr, "done.\n" );
440 for( i=0; i<njob; i++ ) if( kozoarivec[i] ) nkozo++;
444 // if( nlenmax > 30000 )
445 if( nlenmax > 50000 ) // version >= 6.823
450 fprintf( stderr, "\nnlenmax=%d, nagasugi!\n", nlenmax );
455 fprintf( stderr, "\nnevermemsave=1, nlenmax=%d, nagasugi!\n", nlenmax );
459 if( !constraint && !nevermemsave && alg != 'M' )
461 fprintf( stderr, "\nnlenmax=%d, Switching to the memsave mode\n", nlenmax );
467 Read( name, nlen, seq_g );
469 readData_pointer( infp, name, nlen, seq_g );
474 for( i=0; i<njob; i++ )
480 for( i=0; i<njob; i++ )
482 identity *= ( nlen[i] == nlen[0] );
486 fprintf( stderr, "Input pre-aligned data\n" );
489 constants( njob, seq_g );
492 fprintf( stderr, "penalty = %d\n", penalty );
493 fprintf( stderr, "penalty_ex = %d\n", penalty_ex );
494 fprintf( stderr, "offset = %d\n", offset );
504 writePre( njob, name, nlen, seq_g, 1 );
516 c = seqcheck( seq_g );
519 fprintf( stderr, "Illegal character %c\n", c );
522 commongappick( njob, seq_g );
524 if( rnakozo && rnaprediction == 'm' )
526 singlerna = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
527 prep = fopen( "hat4", "r" );
528 if( prep == NULL ) ErrorExit( "Make hat4 using mccaskill." );
529 fprintf( stderr, "Loading 'hat4' ... " );
531 for( i=0; i<njob; i++ )
533 gappick0( nogap1seq[0], seq_g[i] );
534 nogaplen = strlen( nogap1seq[0] );
535 singlerna[i] = (RNApair **)calloc( nogaplen, sizeof( RNApair * ) );
536 for( j=0; j<nogaplen; j++ )
538 singlerna[i][j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
539 singlerna[i][j][0].bestpos = -1;
540 singlerna[i][j][0].bestscore = -1.0;
542 readmccaskill( prep, singlerna[i], nogaplen );
545 fprintf( stderr, "\ndone.\n" );
555 prep = fopen( "hat2", "r" );
556 if( !prep ) ErrorExit( "Make hat2." );
557 readhat2_pointer( prep, njob, name, eff );
560 for( i=0; i<njob-1; i++ )
562 for( j=i+1; j<njob; j++ )
564 printf( " %f", eff[i][j] );
570 veryfastsupg_double_loadtree( njob, eff, topol, len );
571 else if( intop ) // v6.528 deha if( intop ) dattanode intree ga mukou datta.
572 veryfastsupg_double_loadtop( njob, eff, topol, len );
573 else if( treemethod == 'X' || treemethod == 'E' || treemethod == 'q' )
574 // veryfastsupg_double_outtree( njob, eff, topol, len, name );
575 fixed_musclesupg_double_treeout( njob, eff, topol, len, name );
576 else if( treemethod == 'n' )
577 nj( njob, eff, topol, len );
578 else if( treemethod == 's' )
579 spg( njob, eff, topol, len );
580 else if( treemethod == 'p' )
581 upg2( njob, eff, topol, len );
582 else ErrorExit( "Incorrect treemethod.\n" );
585 printf( "utree = %d\n", utree );
589 fprintf( stderr, "\n" );
590 if( ( !utree && kobetsubunkatsu ) || constraint || !bunkatsu )
594 anchors[1] =strlen( seq_g[0] );
599 nseg = searchAnchors( njob, seq_g, segment );
601 fprintf( stderr, "### nseg = %d\n", nseg );
602 fprintf( stderr, "### seq_g[0] = %s\n", seq_g[0] );
603 fprintf( stderr, "### nlenmax = %d\n", nlenmax );
604 fprintf( stderr, "### len = %d\n", strlen( seq_g[0] ) );
609 for( i=0; i<nseg; i++ ) anchors[i+1] = segment[i].center;
610 anchors[nseg+1] = strlen( seq_g[0] );
614 for( i=0; i<nseg; i++ )
615 fprintf( stderr, "anchor[%d] = %d\n", i, anchors[i] );
619 for( i=0; i<njob; i++ ) res_g[i][0] = 0;
621 for( iseg=0; iseg<nseg-1; iseg++ )
623 int tmplen = anchors[iseg+1]-anchors[iseg];
624 int pos = strlen( res_g[0] );
625 for( j=0; j<njob; j++ )
627 strncpy( seq[j], seq_g[j], tmplen );
632 fprintf( stderr, "Segment %3d/%3d %4d-%4d\n", iseg+1, nseg-1, pos+1, pos+1+tmplen );
634 fprintf( trap_g, "Segment %3d/%3d %4d-%4d\n", iseg+1, nseg-1, pos+1, pos+1+tmplen );
637 returnvalue = TreeDependentIteration( njob, name, nlen, seq, bseq, topol, len, alloclen, localhomtable, singlerna, nkozo, kozoarivec );
639 for( i=0; i<njob; i++ )
640 strcat( res_g[i], bseq[i] );
642 FreeCharMtx( seq_g_bk );
644 FreeDoubleMtx( len );
645 FreeDoubleMtx( eff );
646 if( constraint ) FreeLocalHomTable( localhomtable, njob );
649 Write( stdout, njob, name, nlen, bseq );
652 fprintf( stderr, "done\n" );
653 fprintf( trap_g, "done\n" );
658 writePre( njob, name, nlen, res_g, 1 );
660 writeData( stdout, njob, name, nlen, res_g, 1 );
669 signed int main( int argc, char *argv[] )
676 gets( b ); njob = atoi( b );
680 if( strstr( b, "ayhoff" ) ) scoremtx = 1;
681 else if( strstr( b, "dna" ) || strstr( b, "DNA" ) ) scoremtx = -1;
682 else if( strstr( b, "M-Y" ) || strstr( b, "iyata" ) ) scoremtx = 2;
685 if( strstr( b, "constraint" ) ) cnst = 1;
692 if( !strncmp( b, a, 1 ) )
694 gets( b ); nlen[i] = atoi( b );
695 if( nlen[i] > nlenmax ) nlenmax = nlen[i];
699 if( nlenmax > N || njob > M )
701 fprintf( stderr, "ERROR in main\n" );
708 value = main1( nlen, argc, argv );