+++ /dev/null
-#include "mltaln.h"
-
-
-#define TREE 1
-#define PICKSIZE 50 // must be >= 3
-#define WEIGHT 0
-#define TOKYORIPARA 0.70 // 0.70
-#define TOKYORIPARA_A 0.70 // changed
-#define LENFAC 1
-#define HUKINTOTREE 1
-#define DIANA 0
-#define MAX6DIST 10.0
-
-// kouzoutai ni sasareru pointer ha static
-
-#define DEBUG 0
-#define IODEBUG 0
-#define SCOREOUT 0
-
-#define END_OF_VEC -1
-
-static char *fastapath;
-static int doalign;
-static int fromaln;
-static int uselongest;
-static int treeout;
-static int classsize;
-static int picksize;
-static int maxl;
-static int tsize;
-static int reorder;
-static int pid;
-static int maxdepth = 0;
-static double tokyoripara;
-
-static double lenfaca, lenfacb, lenfacc, lenfacd;
-#define PLENFACA 0.01
-#define PLENFACB 10000
-#define PLENFACC 10000
-#define PLENFACD 0.1
-#define DLENFACA 0.01
-#define DLENFACB 2500
-#define DLENFACC 2500
-#define DLENFACD 0.1
-
-static char datafile[1000];
-static char queryfile[1000];
-static char resultfile[1000];
-
-typedef struct _scores
-{
- double score;
- int selfscore;
- int orilen;
- int *pointt;
- int numinseq;
- char *name;
-// char *seq; // reallo
-// char **seqpt;
- int shimon;
-} Scores;
-
-int intcompare( const int *a, const int *b )
-{
- return( *a - *b );
-}
-
-int lcompare( const Scores *a, const Scores *b )
-{
- if( a->orilen < b->orilen ) return 1;
- else if( a->orilen > b->orilen ) return -1;
- else return 0;
-}
-
-int dcompare( const Scores *a, const Scores *b )
-{
- if( a->score > b->score ) return 1;
- else if( a->score < b->score ) return -1;
- else
- {
- if( a->selfscore < b->selfscore ) return 1;
- else if( a->selfscore > b->selfscore ) return -1;
- else
- {
- if( a->orilen < b->orilen ) return 1;
- else if( a->orilen > b->orilen ) return -1;
- else return 0;
- }
- }
-}
-
-#if 0
-static void gappickandx0( char *out, char *in )
-{
- char c;
- if( scoremtx == -1 )
- {
- while( *in )
- {
- if( (c=*in++) == '-' )
- ;
- else if( c == 'u' )
- *out++ = 't';
- else if( amino_n[c] < 4 && amino_n[c] > -1 )
- *out++ = c;
- else
- *out++ = 'n';
- }
- }
- else
- {
- while( *in )
- {
- if( (c=*in++) == '-' )
- ;
- else if( amino_n[c] < 20 && amino_n[c] > -1 )
- *out++ = c;
- else
- *out++ = 'X';
- }
- }
- *out = 0;
-}
-
-static int getkouho( int *pickkouho, double prob, int nin, Scores *scores, char **seq ) // 0 < prob < 1
-{
- int nkouho = 0;
- int i, j;
- int *iptr = pickkouho;
- for( i=1; i<nin; i++ )
- {
- if( ( nkouho==0 || rnd() < prob ) && ( scores[i].shimon != scores->shimon || strcmp( seq[scores->numinseq], seq[scores[i].numinseq] ) ) )
- {
-#if 0
- for( j=0; j<nkouho; j++ )
- {
- if( scores[i].shimon == scores[pickkouho[j]].shimon || !strcmp( seq[scores[pickkouho[j]].numinseq], seq[scores[i].numinseq] ) )
- break;
- }
- if( j == nkouho )
-#endif
- {
- *iptr++ = i;
- nkouho++;
-// fprintf( stderr, "ok! nkouho=%d\n", nkouho );
- }
- }
- else
- {
- ;
-// fprintf( stderr, "no! %d-%d\n", 0, scores[i].numinseq );
- }
- }
- fprintf( stderr, "\ndone\n\n" );
- return nkouho;
-}
-
-#endif
-
-static void getfastascoremtx( int **tmpaminodis )
-{
- FILE *qfp;
- FILE *dfp;
- FILE *rfp;
- int i, j;
- char aa;
- int slen;
- int res;
- char com[10000];
- static char *tmpseq;
- static char *tmpname;
- double *resvec;
-
- if( scoremtx == -1 )
- {
- tmpaminodis['a']['a'] = 5;
- tmpaminodis['g']['g'] = 5;
- tmpaminodis['c']['c'] = 5;
- tmpaminodis['t']['t'] = 5;
- tmpaminodis['n']['n'] = -1;
-
- return;
- }
-
-
- tmpseq = calloc( 2000, sizeof( char ) );
- tmpname = calloc( B, sizeof( char ) );
- resvec = calloc( 1, sizeof( double ) );
-
-// fprintf( stderr, "xformatting .. " );
- dfp = fopen( datafile, "w" );
- if( !dfp ) ErrorExit( "Cannot open datafile." );
- sprintf( tmpname, ">+===========+%d ", 0 );
- strcpy( tmpseq, "AAAAAAXXXXXX" );
- strcat( tmpseq, "CCCCCCXXXXXX" );
- strcat( tmpseq, "DDDDDDXXXXXX" );
- strcat( tmpseq, "EEEEEEXXXXXX" );
- strcat( tmpseq, "FFFFFFXXXXXX" );
- strcat( tmpseq, "GGGGGGXXXXXX" );
- strcat( tmpseq, "HHHHHHXXXXXX" );
- strcat( tmpseq, "IIIIIIXXXXXX" );
- strcat( tmpseq, "KKKKKKXXXXXX" );
- strcat( tmpseq, "LLLLLLXXXXXX" );
- strcat( tmpseq, "MMMMMMXXXXXX" );
- strcat( tmpseq, "NNNNNNXXXXXX" );
- strcat( tmpseq, "PPPPPPXXXXXX" );
- strcat( tmpseq, "QQQQQQXXXXXX" );
- strcat( tmpseq, "RRRRRRXXXXXX" );
- strcat( tmpseq, "SSSSSSXXXXXX" );
- strcat( tmpseq, "TTTTTTXXXXXX" );
- strcat( tmpseq, "VVVVVVXXXXXX" );
- strcat( tmpseq, "WWWWWWXXXXXX" );
- strcat( tmpseq, "YYYYYYXXXXXX" );
- slen = strlen( tmpseq );
- writeData_pointer( dfp, 1, &tmpname, &slen, &tmpseq );
- fclose( dfp );
- fprintf( stderr, "done.\n" );
-
- for( i=0; i<20; i++ )
- {
- aa = amino[i];
-// fprintf( stderr, "checking %c\n", aa );
- *tmpseq = 0;
- sprintf( tmpname, ">+===========+%d ", 0 );
- for( j=0; j<6; j++ )
- sprintf( tmpseq+strlen( tmpseq ), "%c", aa );
- qfp = fopen( queryfile, "w" );
- if( !qfp ) ErrorExit( "Cannot open queryfile." );
- writeData_pointer( qfp, 1, &tmpname, &slen, &tmpseq );
- fclose( qfp );
-
- if( scoremtx == -1 )
- sprintf( com, "%s -z3 -m10 -n -Q -H -b%d -E%d -d%d %s %s %d > %s", fastapath, M, M, 0, queryfile, datafile, 6, resultfile );
- else
- sprintf( com, "%s -z3 -m10 -p -Q -H -b%d -E%d -d%d %s %s %d > %s", fastapath, M, M, 0, queryfile, datafile, 2, resultfile );
- res = system( com );
- if( res )
- {
- fprintf( stderr, "error in %s", fastapath );
- exit( 1 );
- }
-
- rfp = fopen( resultfile, "r" );
- if( rfp == NULL )
- ErrorExit( "file 'fasta.$$' does not exist\n" );
- res = ReadFasta34m10_scoreonly( rfp, resvec, 1 );
- fprintf( stderr, "%c: %f\n", 'A'+i, *resvec/6 );
- fclose( rfp );
- if( ( (int)*resvec % 6 ) > 0.0 )
- {
- fprintf( stderr, "Error in blast, *resvec=%f\n", *resvec );
- fprintf( stderr, "Error in blast, *resvec/6=%f\n", *resvec/6 );
- exit( 1 );
- }
- tmpaminodis[(int)aa][(int)aa] = (int)( *resvec / 6 );
-// fprintf( stderr, "*resvec=%f, tmpaminodis[aa][aa] = %d\n", *resvec, tmpaminodis[aa][aa] );
- }
- tmpaminodis['X']['X'] = -1;
- free( tmpname );
- free( tmpseq );
- free( resvec );
-}
-
-#if 0
-static void getblastscoremtx( int **tmpaminodis )
-{
- FILE *qfp;
- FILE *dfp;
- FILE *rfp;
- int i, j;
- char aa;
- int slen;
- int res;
- char com[10000];
- static char *tmpseq;
- static char *tmpname;
- double *resvec;
-
- if( scoremtx == -1 )
- {
- tmpaminodis['a']['a'] = 1;
- tmpaminodis['g']['g'] = 1;
- tmpaminodis['c']['c'] = 1;
- tmpaminodis['t']['t'] = 1;
-
- return;
- }
-
-
- tmpseq = calloc( 2000, sizeof( char ) );
- tmpname = calloc( B, sizeof( char ) );
- resvec = calloc( 1, sizeof( double ) );
-
-// fprintf( stderr, "xformatting .. " );
- dfp = fopen( datafile, "w" );
- if( !dfp ) ErrorExit( "Cannot open datafile." );
- sprintf( tmpname, "\0", i ); // BUG!!
- strcpy( tmpseq, "AAAAAAXXXXXX" );
- strcat( tmpseq, "CCCCCCXXXXXX" );
- strcat( tmpseq, "DDDDDDXXXXXX" );
- strcat( tmpseq, "EEEEEEXXXXXX" );
- strcat( tmpseq, "FFFFFFXXXXXX" );
- strcat( tmpseq, "GGGGGGXXXXXX" );
- strcat( tmpseq, "HHHHHHXXXXXX" );
- strcat( tmpseq, "IIIIIIXXXXXX" );
- strcat( tmpseq, "KKKKKKXXXXXX" );
- strcat( tmpseq, "LLLLLLXXXXXX" );
- strcat( tmpseq, "MMMMMMXXXXXX" );
- strcat( tmpseq, "NNNNNNXXXXXX" );
- strcat( tmpseq, "PPPPPPXXXXXX" );
- strcat( tmpseq, "QQQQQQXXXXXX" );
- strcat( tmpseq, "RRRRRRXXXXXX" );
- strcat( tmpseq, "SSSSSSXXXXXX" );
- strcat( tmpseq, "TTTTTTXXXXXX" );
- strcat( tmpseq, "VVVVVVXXXXXX" );
- strcat( tmpseq, "WWWWWWXXXXXX" );
- strcat( tmpseq, "YYYYYYXXXXXX" );
- slen = strlen( tmpseq );
- writeData_pointer( dfp, 1, &tmpname, &slen, &tmpseq );
- fclose( dfp );
- if( scoremtx == -1 )
- sprintf( com, "formatdb -p f -i %s -o F", datafile );
- else
- sprintf( com, "formatdb -i %s -o F", datafile );
- system( com );
- fprintf( stderr, "done.\n" );
-
- for( i=0; i<20; i++ )
- {
- aa = amino[i];
- fprintf( stderr, "checking %c\n", aa );
- *tmpseq = 0;
- for( j=0; j<6; j++ )
- sprintf( tmpseq+strlen( tmpseq ), "%c", aa );
- qfp = fopen( queryfile, "w" );
- if( !qfp ) ErrorExit( "Cannot open queryfile." );
- writeData_pointer( qfp, 1, &tmpname, &slen, &tmpseq );
- fclose( qfp );
-
- sprintf( com, "blastall -b %d -G 10 -E 1 -e 1e10 -p blastp -m 7 -i %s -d %s > %s\0", 1, queryfile, datafile, resultfile );
- res = system( com );
- if( res )
- {
- fprintf( stderr, "error in %s", "blastall" );
- exit( 1 );
- }
-
- rfp = fopen( resultfile, "r" );
- if( rfp == NULL )
- ErrorExit( "file 'fasta.$$' does not exist\n" );
- res = ReadBlastm7_scoreonly( rfp, resvec, 1 );
- fprintf( stdout, "%c: %f\n", 'A'+i, *resvec/6 );
- fclose( rfp );
- if( ( (int)*resvec % 6 ) > 0.0 )
- {
- fprintf( stderr, "Error in blast, *resvec=%f\n", *resvec );
- fprintf( stderr, "Error in blast, *resvec/6=%f\n", *resvec/6 );
- exit( 1 );
- }
- tmpaminodis[aa][aa] = (int)( *resvec / 6 );
- }
- tmpaminodis['X']['X'] = 0;
- free( tmpname );
- free( tmpseq );
- free( resvec );
-
-}
-#endif
-
-static double *callfasta( char **seq, Scores *scores, int nin, int *picks, int query, int rewritedata )
-{
- double *val;
- FILE *qfp;
- FILE *dfp;
- FILE *rfp;
- int i;
- char com[10000];
- static char datafile[1000];
- static char queryfile[1000];
- static char resultfile[1000];
- static int pid;
- static char *tmpseq;
- static char *tmpname;
- int slen;
- int res;
- static Scores *scoresbk = NULL;
- static int ninbk = 0;
-
- if( pid == 0 )
- {
- pid = (int)getpid();
- sprintf( datafile, "/tmp/data-%d", pid );
- sprintf( queryfile, "/tmp/query-%d", pid );
- sprintf( resultfile, "/tmp/fasta-%d", pid );
-
- tmpseq = calloc( nlenmax+1, sizeof( char ) );
- tmpname = calloc( B+1, sizeof( char ) );
- }
-
- val = calloc( nin, sizeof( double ) );
-// fprintf( stderr, "nin=%d, q=%d\n", nin, query );
-
- if( rewritedata )
- {
- scoresbk = scores;
- ninbk = nin;
-// fprintf( stderr, "\nformatting .. " );
- dfp = fopen( datafile, "w" );
- if( !dfp ) ErrorExit( "Cannot open datafile." );
- if( picks == NULL ) for( i=0; i<nin; i++ )
- {
-// fprintf( stderr, "i=%d / %d / %d\n", i, nin, njob );
-// fprintf( stderr, "nlenmax = %d\n", nlenmax );
-// fprintf( stderr, "scores[i].orilen = %d\n", scores[i].orilen );
-// fprintf( stderr, "strlen( seq[scores[i].numinseq] = %d\n", strlen( seq[scores[i].numinseq] ) );
- gappick0( tmpseq, seq[scores[i].numinseq] );
- sprintf( tmpname, ">+===========+%d ", i );
- slen = scores[i].orilen;
- writeData_pointer( dfp, 1, &tmpname, &slen, &tmpseq );
- }
- else for( i=0; i<nin; i++ )
- {
- gappick0( tmpseq, seq[scores[picks[i]].numinseq] );
- sprintf( tmpname, ">+===========+%d ", i );
- slen = scores[picks[i]].orilen;
- writeData_pointer( dfp, 1, &tmpname, &slen, &tmpseq );
- }
- fclose( dfp );
- }
-
-
- gappick0( tmpseq, seq[scores[query].numinseq] );
- sprintf( tmpname, ">+==========+%d ", 0 );
- slen = scores[query].orilen;
- qfp = fopen( queryfile, "w" );
- if( !qfp ) ErrorExit( "Cannot open queryfile." );
- writeData_pointer( qfp, 1, &tmpname, &slen, &tmpseq );
- fclose( qfp );
-
-// fprintf( stderr, "calling fasta, nin=%d\n", nin );
-
- if( scoremtx == -1 )
- sprintf( com, "%s -z3 -m10 -n -Q -H -b%d -E%d -d%d %s %s %d > %s", fastapath, nin, nin, 0, queryfile, datafile, 6, resultfile );
- else
- sprintf( com, "%s -z3 -m10 -p -Q -H -b%d -E%d -d%d %s %s %d > %s", fastapath, nin, nin, 0, queryfile, datafile, 2, resultfile );
- res = system( com );
- if( res )
- {
- fprintf( stderr, "error in %s", fastapath );
- exit( 1 );
- }
-// fprintf( stderr, "fasta done\n" );
-
-//exit( 1 );
-
- rfp = fopen( resultfile, "r" );
- if( rfp == NULL )
- ErrorExit( "file 'fasta.$$' does not exist\n" );
-
-// fprintf( stderr, "reading fasta\n" );
- if( scoremtx == -1 )
- res = ReadFasta34m10_scoreonly_nuc( rfp, val, nin );
- else
- res = ReadFasta34m10_scoreonly( rfp, val, nin );
-// fprintf( stderr, "done. val[0] = %f\n", val[0] );
-
-
- fclose( rfp );
-
-#if 0
- for( i=0; i<nin; i++ )
- fprintf( stderr, "r[%d-%d] = %f\n", 0, i, val[i] );
- exit( 1 );
-#endif
-
- return( val );
-}
-#if 0
-static double *callblast( char **seq, Scores *scores, int nin, int query, int rewritedata )
-{
- double *val;
- FILE *qfp;
- FILE *dfp;
- FILE *rfp;
- int i, j;
- char com[10000];
- static char datafile[1000];
- static char queryfile[1000];
- static char resultfile[1000];
- static int pid;
- static char *tmpseq;
- static char *tmpname;
- char *seqptr;
- int slen;
- int res;
- static Scores *scoresbk = NULL;
- static int ninbk = 0;
-
- if( pid == 0 )
- {
- pid = (int)getpid();
- sprintf( datafile, "/tmp/data-%d\0", pid );
- sprintf( queryfile, "/tmp/query-%d\0", pid );
- sprintf( resultfile, "/tmp/fasta-%d\0", pid );
-
- tmpseq = calloc( nlenmax+1, sizeof( char ) );
- tmpname = calloc( B+1, sizeof( char ) );
- }
-
- val = calloc( nin, sizeof( double ) );
-// fprintf( stderr, "nin=%d, q=%d\n", nin, query );
-
- if( rewritedata )
- {
- scoresbk = scores;
- ninbk = nin;
- fprintf( stderr, "\nformatting .. " );
- dfp = fopen( datafile, "w" );
- if( !dfp ) ErrorExit( "Cannot open datafile." );
- for( i=0; i<nin; i++ )
- {
-// fprintf( stderr, "i=%d / %d / %d\n", i, nin, njob );
-// fprintf( stderr, "nlenmax = %d\n", nlenmax );
-// fprintf( stderr, "scores[i].orilen = %d\n", scores[i].orilen );
-// fprintf( stderr, "strlen( seq[scores[i].numinseq] = %d\n", strlen( seq[scores[i].numinseq] ) );
- gappick0( tmpseq, seq[scores[i].numinseq] );
- sprintf( tmpname, "+===========+%d \0", i );
- slen = scores[i].orilen;
- writeData_pointer( dfp, 1, &tmpname, &slen, &tmpseq );
- }
- fclose( dfp );
-
- if( scoremtx == -1 )
- sprintf( com, "formatdb -p f -i %s -o F", datafile );
- else
- sprintf( com, "formatdb -i %s -o F", datafile );
- system( com );
-// fprintf( stderr, "done.\n" );
- }
-
-
- gappick0( tmpseq, seq[scores[query].numinseq] );
- sprintf( tmpname, "+==========+%d \0", 0 );
- slen = scores[query].orilen;
- qfp = fopen( queryfile, "w" );
- if( !qfp ) ErrorExit( "Cannot open queryfile." );
- writeData_pointer( qfp, 1, &tmpname, &slen, &tmpseq );
- fclose( qfp );
-// fprintf( stderr, "q=%s\n", tmpseq );
-
- fprintf( stderr, "\ncalling blast .. \n" );
- if( scoremtx == -1 )
- sprintf( com, "blastall -b %d -e 1e10 -p blastn -m 7 -i %s -d %s > %s\0", nin, queryfile, datafile, resultfile );
- else
- sprintf( com, "blastall -b %d -G 10 -E 1 -e 1e10 -p blastp -m 7 -i %s -d %s > %s\0", nin, queryfile, datafile, resultfile );
- res = system( com );
- if( res ) ErrorExit( "error in blast" );
-
- rfp = fopen( resultfile, "r" );
- if( rfp == NULL )
- ErrorExit( "file 'fasta.$$' does not exist\n" );
- res = ReadBlastm7_scoreonly( rfp, val, nin );
- fclose( rfp );
-
-#if 0
- for( i=0; i<nin; i++ )
- fprintf( stderr, "r[%d-%d] = %f\n", 0, i, val[i] );
-#endif
-
- return( val );
-}
-#endif
-
-#if 0
-static void selhead( int *ar, int n )
-{
- int min = *ar;
- int *minptr = ar;
- int *ptr = ar;
- int tmp;
- n--;
- ar++;
- while( n-- )
- {
- if( ( tmp = *ptr++ ) < min )
- {
- min = tmp;
- minptr = ptr;
- }
- }
- if( minptr != ar )
- {
- tmp = *ar;
- *ar = min;
- *minptr = tmp;
- }
- return;
-}
-#endif
-
-void arguments( int argc, char *argv[] )
-{
- int c;
-
- doalign = 0;
- fromaln = 0;
- treeout = 0;
- uselongest = 1;
- reorder = 1;
- nevermemsave = 0;
- inputfile = NULL;
- fftkeika = 0;
- constraint = 0;
- nblosum = 62;
- fmodel = 0;
- calledByXced = 0;
- devide = 0;
- use_fft = 0;
- force_fft = 0;
- fftscore = 1;
- fftRepeatStop = 0;
- fftNoAnchStop = 0;
- weight = 3;
- utree = 1;
- tbutree = 1;
- refine = 0;
- check = 1;
- cut = 0.0;
- disp = 0;
- outgap = 1;
- alg = 'A';
- mix = 0;
- tbitr = 0;
- scmtd = 5;
- tbweight = 0;
- tbrweight = 3;
- checkC = 0;
- treemethod = 'X';
- contin = 0;
- scoremtx = 1;
- kobetsubunkatsu = 0;
- dorp = NOTSPECIFIED;
- ppenalty = -1530;
- ppenalty_ex = NOTSPECIFIED;
- poffset = -123;
- kimuraR = NOTSPECIFIED;
- pamN = NOTSPECIFIED;
- geta2 = GETA2;
- fftWinSize = NOTSPECIFIED;
- fftThreshold = NOTSPECIFIED;
- TMorJTT = JTT;
- classsize = NOTSPECIFIED;
- picksize = NOTSPECIFIED;
- tokyoripara = NOTSPECIFIED;
-
- while( --argc > 0 && (*++argv)[0] == '-' )
- {
- while ( ( c = *++argv[0] ) )
- {
- switch( c )
- {
- case 'p':
- picksize = atoi( *++argv );
- fprintf( stderr, "picksize = %d\n", picksize );
- --argc;
- goto nextoption;
- case 's':
- classsize = atoi( *++argv );
- fprintf( stderr, "groupsize = %d\n", classsize );
- --argc;
- goto nextoption;
- case 'i':
- inputfile = *++argv;
- fprintf( stderr, "inputfile = %s\n", inputfile );
- --argc;
- goto nextoption;
- case 'f':
- ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
-// fprintf( stderr, "ppenalty = %d\n", ppenalty );
- --argc;
- goto nextoption;
- case 'g':
- ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
- fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
- --argc;
- goto nextoption;
- case 'h':
- poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
-// fprintf( stderr, "poffset = %d\n", poffset );
- --argc;
- goto nextoption;
- case 'k':
- kimuraR = atoi( *++argv );
- fprintf( stderr, "kimuraR = %d\n", kimuraR );
- --argc;
- goto nextoption;
- case 'b':
- nblosum = atoi( *++argv );
- scoremtx = 1;
-// fprintf( stderr, "blosum %d\n", nblosum );
- --argc;
- goto nextoption;
- case 'j':
- pamN = atoi( *++argv );
- scoremtx = 0;
- TMorJTT = JTT;
- fprintf( stderr, "jtt %d\n", pamN );
- --argc;
- goto nextoption;
- case 'm':
- pamN = atoi( *++argv );
- scoremtx = 0;
- TMorJTT = TM;
- fprintf( stderr, "tm %d\n", pamN );
- --argc;
- goto nextoption;
- case 'T':
- tokyoripara = (double)atof( *++argv );
- --argc;
- goto nextoption;
- case 'l':
- uselongest = 0;
- break;
-#if 1
- case 'a':
- fmodel = 1;
- break;
-#endif
- case 'S':
- doalign = 'f';
- break;
- case 'Z':
- fromaln = 1;
- break;
- case 'L':
- doalign = 1;
- break;
- case 'x':
- reorder = 0;
- break;
- case 't':
- treeout = 1;
- break;
- case 'r':
- fmodel = -1;
- break;
- case 'D':
- dorp = 'd';
- break;
- case 'P':
- dorp = 'p';
- break;
- case 'e':
- fftscore = 0;
- break;
- case 'O':
- fftNoAnchStop = 1;
- break;
-#if 0
- case 'R':
- fftRepeatStop = 1;
- break;
- case 'Q':
- calledByXced = 1;
- break;
- case 'a':
- alg = 'a';
- break;
-#endif
- case 'R':
- alg = 'R';
- break;
- case 'Q':
- alg = 'Q';
- break;
- case 'A':
- alg = 'A';
- break;
- case 'N':
- nevermemsave = 1;
- break;
- case 'M':
- alg = 'M';
- break;
- case 'C':
- alg = 'C';
- break;
- case 'F':
- use_fft = 1;
- break;
- case 'G':
- use_fft = 1;
- force_fft = 1;
- break;
- case 'v':
- tbrweight = 3;
- break;
- case 'd':
- disp = 1;
- break;
- case 'o':
- outgap = 0;
- break;
- case 'J':
- tbutree = 0;
- break;
- case 'X':
- treemethod = 'X'; // mix
- break;
- case 'E':
- treemethod = 'E'; // upg (average)
- break;
- case 'q':
- treemethod = 'q'; // minimum
- break;
- case 'z':
- fftThreshold = atoi( *++argv );
- --argc;
- goto nextoption;
- case 'w':
- fftWinSize = atoi( *++argv );
- --argc;
- goto nextoption;
- default:
- fprintf( stderr, "illegal option %c\n", c );
- argc = 0;
- break;
- }
- }
- nextoption:
- ;
- }
- if( argc == 1 )
- {
- cut = atof( (*argv) );
- argc--;
- }
- if( argc != 0 )
- {
- fprintf( stderr, "options: Check source file !\n" );
- exit( 1 );
- }
- if( tbitr == 1 && outgap == 0 )
- {
- fprintf( stderr, "conflicting options : o, m or u\n" );
- exit( 1 );
- }
- if( alg == 'C' && outgap == 0 )
- {
- fprintf( stderr, "conflicting options : C, o\n" );
- exit( 1 );
- }
-}
-
-static int maxl;
-static int tsize;
-
-int seq_grp_nuc( int *grp, char *seq )
-{
- int tmp;
- int *grpbk = grp;
- while( *seq )
- {
- tmp = amino_grp[(int)*seq++];
- if( tmp < 4 )
- *grp++ = tmp;
- else
- fprintf( stderr, "WARNING : Unknown character %c\r", *(seq-1) );
- }
- *grp = END_OF_VEC;
- return( grp-grpbk );
-}
-
-int seq_grp( int *grp, char *seq )
-{
- int tmp;
- int *grpbk = grp;
- while( *seq )
- {
- tmp = amino_grp[(int)*seq++];
- if( tmp < 6 )
- *grp++ = tmp;
- else
- fprintf( stderr, "WARNING : Unknown character %c\r", *(seq-1) );
- }
- *grp = END_OF_VEC;
- return( grp-grpbk );
-}
-
-void makecompositiontable_p( short *table, int *pointt )
-{
- int point;
-
- while( ( point = *pointt++ ) != END_OF_VEC )
- table[point]++;
-}
-
-int commonsextet_p( short *table, int *pointt )
-{
- int value = 0;
- short tmp;
- int point;
- static short *memo = NULL;
- static int *ct = NULL;
- static int *cp;
-
- if( !memo )
- {
- memo = (short *)calloc( tsize, sizeof( short ) );
- if( !memo ) ErrorExit( "Cannot allocate memo\n" );
- ct = (int *)calloc( MIN( maxl, tsize )+1, sizeof( int ) );
- if( !ct ) ErrorExit( "Cannot allocate memo\n" );
- }
-
- cp = ct;
- while( ( point = *pointt++ ) != END_OF_VEC )
- {
- tmp = memo[point]++;
- if( tmp < table[point] )
- value++;
- if( tmp == 0 ) *cp++ = point;
- }
- *cp = END_OF_VEC;
-
- cp = ct;
- while( *cp != END_OF_VEC )
- memo[*cp++] = 0;
-
- return( value );
-}
-
-void makepointtable_nuc( int *pointt, int *n )
-{
- int point;
- register int *p;
-
- p = n;
- point = *n++ * 1024;
- point += *n++ * 256;
- point += *n++ * 64;
- point += *n++ * 16;
- point += *n++ * 4;
- point += *n++;
- *pointt++ = point;
-
- while( *n != END_OF_VEC )
- {
- point -= *p++ * 1024;
- point *= 4;
- point += *n++;
- *pointt++ = point;
- }
- *pointt = END_OF_VEC;
-}
-
-void makepointtable( int *pointt, int *n )
-{
- int point;
- register int *p;
-
- p = n;
- point = *n++ * 7776;
- point += *n++ * 1296;
- point += *n++ * 216;
- point += *n++ * 36;
- point += *n++ * 6;
- point += *n++;
- *pointt++ = point;
-
- while( *n != END_OF_VEC )
- {
- point -= *p++ * 7776;
- point *= 6;
- point += *n++;
- *pointt++ = point;
- }
- *pointt = END_OF_VEC;
-}
-
-#if 1
-static void pairalign( int nseq, int *nlen, char **seq, int *mem1, int *mem2, double *weight, int *alloclen )
-{
- int l, len1, len2;
- int clus1, clus2;
- float pscore, tscore;
- static int *fftlog;
- static char *indication1, *indication2;
- static double *effarr1 = NULL;
- static double *effarr2 = NULL;
- static char **mseq1, **mseq2;
- float dumfl = 0.0;
- int ffttry;
- int m1, m2;
-#if 0
- int i, j;
-#endif
-
-
- if( effarr1 == NULL )
- {
- fftlog = AllocateIntVec( nseq );
- effarr1 = AllocateDoubleVec( nseq );
- effarr2 = AllocateDoubleVec( nseq );
- indication1 = AllocateCharVec( 150 );
- indication2 = AllocateCharVec( 150 );
- mseq1 = AllocateCharMtx( nseq, 0 );
- mseq2 = AllocateCharMtx( nseq, 0 );
- for( l=0; l<nseq; l++ ) fftlog[l] = 1;
- }
-
- tscore = 0.0;
- m1 = mem1[0];
- m2 = mem2[0];
- len1 = strlen( seq[m1] );
- len2 = strlen( seq[m2] );
- if( *alloclen < len1 + len2 )
- {
- fprintf( stderr, "\nReallocating.." );
- *alloclen = ( len1 + len2 ) + 1000;
- ReallocateCharMtx( seq, nseq, *alloclen + 10 );
- fprintf( stderr, "done. *alloclen = %d\n", *alloclen );
- }
-
-#if WEIGHT
- clus1 = fastconjuction_noname( mem1, seq, mseq1, effarr1, weight, indication1 );
- clus2 = fastconjuction_noname( mem2, seq, mseq2, effarr2, weight, indication2 );
-#else
- clus1 = fastconjuction_noweight( mem1, seq, mseq1, effarr1, indication1 );
- clus2 = fastconjuction_noweight( mem2, seq, mseq2, effarr2, indication2 );
-#endif
-
-#if 0
- for( i=0; i<clus1; i++ )
- fprintf( stderr, "in p seq[%d] = %s\n", mem1[i], seq[mem1[i]] );
- for( i=0; i<clus2; i++ )
- fprintf( stderr, "in p seq[%d] = %s\n", mem2[i], seq[mem2[i]] );
-#endif
-
-#if 0
- fprintf( stderr, "group1 = %.66s", indication1 );
- if( strlen( indication1 ) > 66 ) fprintf( stderr, "..." );
- fprintf( stderr, "\n" );
- fprintf( stderr, "group2 = %.66s", indication2 );
- if( strlen( indication2 ) > 66 ) fprintf( stderr, "..." );
- fprintf( stderr, "\n" );
-#endif
-
-// fprintf( stdout, "mseq1 = %s\n", mseq1[0] );
-// fprintf( stdout, "mseq2 = %s\n", mseq2[0] );
-
- if( !nevermemsave && ( alg != 'M' && ( len1 > 10000 || len2 > 10000 ) ) )
- {
- fprintf( stderr, "\nlen1=%d, len2=%d, Switching to the memsave mode\n", len1, len2 );
- alg = 'M';
- if( commonIP ) FreeIntMtx( commonIP );
- commonAlloc1 = 0;
- commonAlloc2 = 0;
- }
-
- if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 );
- else ffttry = 0;
-
- if( force_fft || ( use_fft && ffttry ) )
- {
- fprintf( stderr, "\bf" );
- if( alg == 'M' )
- {
- fprintf( stderr, "\bm" );
-// fprintf( stderr, "%d-%d", clus1, clus2 );
- pscore = Falign_udpari_long( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1 );
- }
- else
- {
-// fprintf( stderr, "%d-%d", clus1, clus2 );
- pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1, NULL, 0, NULL );
- }
- }
- else
- {
- fprintf( stderr, "\bd" );
- fftlog[m1] = 0;
- switch( alg )
- {
- case( 'a' ):
- pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen );
- break;
- case( 'M' ):
- fprintf( stderr, "\bm" );
-// fprintf( stderr, "%d-%d", clus1, clus2 );
- pscore = MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
- break;
- case( 'Q' ):
- if( clus1 == 1 && clus2 == 1 )
- {
-// fprintf( stderr, "%d-%d", clus1, clus2 );
- pscore = G__align11( mseq1, mseq2, *alloclen, outgap, outgap );
- }
- else
- {
-// fprintf( stderr, "%d-%d", clus1, clus2 );
- pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
- }
- break;
- case( 'A' ):
- if( clus1 == 1 && clus2 == 1 )
- {
-// fprintf( stderr, "%d-%d", clus1, clus2 );
- pscore = G__align11( mseq1, mseq2, *alloclen, outgap, outgap );
- }
- else
- {
-// fprintf( stderr, "%d-%d", clus1, clus2 );
- pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
- }
- break;
- default:
- ErrorExit( "ERROR IN SOURCE FILE" );
- }
- }
-#if SCOREOUT
- fprintf( stderr, "score = %10.2f\n", pscore );
-#endif
- nlen[m1] = 0.5 * ( nlen[m1] + nlen[m2] );
- return;
-}
-#endif
-
-#if 0
-static void treebase( int nseq, int *nlen, char **aseq, double *eff, int nalign, int ***topol, int *alloclen ) // topol
-{
- int i, l;
- int nlim;
- int clus1, clus2;
-
- nlim = nalign-1;
- for( l=0; l<nlim; l++ )
- {
- fprintf( stderr, "in treebase, l = %d\n", l );
- fprintf( stderr, "aseq[0] = %s\n", aseq[0] );
- fprintf( stderr, "aseq[topol[l][0][0]] = %s\n", aseq[topol[l][0][0]] );
- pairalign( nseq, nlen, aseq, topol[l][0], topol[l][1], eff, alloclen );
- free( topol[l][0] );
- free( topol[l][1] );
- }
-}
-#endif
-
-static void WriteOptions( FILE *fp )
-{
-
- if( dorp == 'd' ) fprintf( fp, "DNA\n" );
- else
- {
- if ( scoremtx == 0 ) fprintf( fp, "JTT %dPAM\n", pamN );
- else if( scoremtx == 1 ) fprintf( fp, "BLOSUM %d\n", nblosum );
- else if( scoremtx == 2 ) fprintf( fp, "M-Y\n" );
- }
- fprintf( stderr, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
- if( use_fft ) fprintf( fp, "FFT on\n" );
-
- fprintf( fp, "tree-base method\n" );
- if( tbrweight == 0 ) fprintf( fp, "unweighted\n" );
- else if( tbrweight == 3 ) fprintf( fp, "clustalw-like weighting\n" );
- if( tbitr || tbweight )
- {
- fprintf( fp, "iterate at each step\n" );
- if( tbitr && tbrweight == 0 ) fprintf( fp, " unweighted\n" );
- if( tbitr && tbrweight == 3 ) fprintf( fp, " reversely weighted\n" );
- if( tbweight ) fprintf( fp, " weighted\n" );
- fprintf( fp, "\n" );
- }
-
- fprintf( fp, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
-
- if( alg == 'a' )
- fprintf( fp, "Algorithm A\n" );
- else if( alg == 'A' )
- fprintf( fp, "Algorithm A+\n" );
- else if( alg == 'S' )
- fprintf( fp, "Apgorithm S\n" );
- else if( alg == 'C' )
- fprintf( fp, "Apgorithm A+/C\n" );
- else
- fprintf( fp, "Unknown algorithm\n" );
-
- if( treemethod == 'X' )
- fprintf( fp, "Tree = UPGMA (mix).\n" );
- else if( treemethod == 'E' )
- fprintf( fp, "Tree = .UPGMA (average)\n" );
- else if( treemethod == 'q' )
- fprintf( fp, "Tree = Minimum linkage.\n" );
- else
- fprintf( fp, "Unknown tree.\n" );
-
- if( use_fft )
- {
- fprintf( fp, "FFT on\n" );
- if( dorp == 'd' )
- fprintf( fp, "Basis : 4 nucleotides\n" );
- else
- {
- if( fftscore )
- fprintf( fp, "Basis : Polarity and Volume\n" );
- else
- fprintf( fp, "Basis : 20 amino acids\n" );
- }
- fprintf( fp, "Threshold of anchors = %d%%\n", fftThreshold );
- fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
- }
- else
- fprintf( fp, "FFT off\n" );
- fflush( fp );
-}
-
-#if 1
-static int splitseq_mq( Scores *scores, int nin, int *nlen, char **seq, char **orialn, char **name, char *inputfile, int uniform, char **tree, int *alloclen, int *order, int *whichgroup, double *weight, int *depthpt, int qinoya )
-{
- int val;
- int ii, jj;
- int *mptr;
- int treelen = 0; // by Mathog, a guess
- static int groupid = 0;
- static int branchid = 0;
- int i, j;
- int selfscore0;
- double **dfromc;
- double **dfromcp;
- float **pickmtx;
- float **yukomtx;
- static short *table1;
- Scores **outs, *ptr;
- int *numin;
- int *tsukau;
- int belongto;
- char **children = NULL; // by Mathog, a guess
- char *tmptree;
- static int *orderpos = NULL;
- int rn;
- int npick;
- int nyuko;
- int *picks;
- int *yukos;
- int *s_p_map;
- int *p_o_map;
- int *s_y_map;
- int *y_o_map;
- int *closeh;
- int nkouho;
- int *pickkouho;
- int *iptr;
- int *jptr;
- int aligned;
- int ***topol;
- int *treeorder;
- int picktmp;
- float **len;
- double minscore;
-// double *minscoreinpick;
- float *hanni;
- double lenfac;
- double longer;
- double shorter;
- static char **mseq1 = NULL;
- static char **mseq2 = NULL;
- double *blastresults = NULL; // by Mathog, a guess
- static int palloclen = 0;
- float maxdist;
-
- if( orderpos == NULL )
- orderpos = order;
- if( palloclen == 0 )
- palloclen = *alloclen * 2;
- if( mseq1 == NULL && doalign == 1 )
- {
- mseq1 = AllocateCharMtx( 1, palloclen );
- mseq2 = AllocateCharMtx( 1, palloclen );
- }
-
- if( nin == 0 )
- {
-#if TREE
- if( treeout )
- {
- *tree = (char *)calloc( 1, sizeof( char ) );
- **tree = 0;
- }
-#endif
- return 1;
- }
-
- if( nin < 2 || uniform == -1 ) // kokodato muda deha nai ga
- {
- fprintf( stderr, "\nLeaf %d / %d ", ++branchid, njob );
-#if 0
- outputfile = AllocateCharVec( strlen( inputfile ) + 100 );
- sprintf( outputfile, "%s-%d", inputfile, branchid );
- if( uniform > 0 )
-// sprintf( outputfile, "%su%d", outputfile, uniform );
- sprintf( outputfile + strlen(outputfile), "u%d", uniform );
- fprintf( stderr, "GROUP %d: %d member(s) (%d) %s\n", branchid, nin, scores[0].numinseq, outputfile );
- outfp = fopen( outputfile, "w" );
- free( outputfile );
- if( outfp == NULL )
- {
- fprintf( stderr, "Cannot open %s\n", outputfile );
- exit( 1 );
- }
- for( j=0; j<nin; j++ )
- fprintf( outfp, ">G%d %s\n%s\n", branchid, scores[j].name+1, seq[scores[j].numinseq] );
- fclose( outfp );
-#endif
-
-
-#if TREE
- if( treeout )
- {
- treelen = 0;
- tmptree = calloc( 100, sizeof( char ) );
- for( j=0; j<nin; j++ )
- {
- treelen += sprintf( tmptree, "%d", scores[j].numinseq+1 );
- }
- free( tmptree );
-
- *tree = (char *)calloc( treelen + nin + 5, sizeof( char ) );
- if( nin > 1 ) **tree = '(';
- else **tree = '\0';
-// **tree = '\0';
- for( j=0; j<nin-1; j++ )
- {
- sprintf( *tree+strlen( *tree ), "%d,", scores[j].numinseq+1 );
- }
- sprintf( *tree+strlen( *tree ), "%d", scores[j].numinseq+1 );
- if( nin > 1 ) strcat( *tree, ")" );
-// fprintf( stdout, "*tree = %s\n", *tree );
- }
-
-#endif
- for( j=0; j<nin; j++ )
- {
- *orderpos++ = scores[j].numinseq;
-// fprintf( stderr, "*order = %d\n", scores[j].numinseq );
- }
-
- return 1;
- }
-
-
-
- if( uselongest )
- {
- i = nin;
- ptr = scores;
- selfscore0 = scores->selfscore;
- belongto = 0;
- while( i-- )
- {
-// fprintf( stderr, "ptr-scores=%d, numinseq = %d, score = %f\n", ptr-scores, ptr->numinseq+1, ptr->score );
- if( ptr->selfscore > selfscore0 )
- {
- selfscore0 = ptr->selfscore;
- belongto = ptr-scores;
- }
- ptr++;
- }
-#if 1
- if( belongto != 0 )
- {
-// fprintf( stderr, "swap %d %s\n<->\n%d %s\n", 0, scores->name, belongto, (scores+belongto)->name );
- ptr = calloc( 1, sizeof( Scores ) );
- *ptr = scores[belongto];
- scores[belongto] = *scores;
- *scores = *ptr;
- free( ptr );
- }
-#endif
- }
- else
- {
- qsort( scores, nin, sizeof( Scores ), (int (*)())lcompare );
- belongto = (int)( 0.5 * nin );
-// fprintf( stderr, "lengths = %d, %d, %d\n", scores->orilen, scores[belongto].orilen, scores[nin-1].orilen );
- if( belongto != 0 )
- {
-// fprintf( stderr, "swap %d %s\n<->\n%d %s\n", 0, scores->name, belongto, (scores+belongto)->name );
- ptr = calloc( 1, sizeof( Scores ) );
- *ptr = scores[belongto];
- scores[belongto] = *scores;
- *scores = *ptr;
- free( ptr );
- }
- }
-
- if( qinoya != scores->numinseq )
-// if( 1 || qinoya != scores->numinseq )
- {
-// fprintf( stdout, "### scores->numinseq = %d, qinoya=%d, depth=%d\n", scores->numinseq, qinoya, *depthpt );
-
-
- if( doalign )
- {
- if( doalign == 'f' )
- {
- blastresults = callfasta( seq, scores, nin, NULL, 0, 1 );
- if( scores->selfscore != (int)blastresults[0] )
- {
- fprintf( stderr, "\n\nWARNING1: selfscore\n" );
- fprintf( stderr, "scores->numinseq = %d\n", scores->numinseq+1 );
- fprintf( stderr, "scores->orilen = %d\n", scores->orilen );
- fprintf( stderr, "scores->selfscore = %d, but blastresults[0] = %f\n", scores->selfscore, blastresults[0] );
-// if( abs( scores->selfscore - (int)blastresults[0] ) > 2 )
-// exit( 1 );
-// scores->selfscore = (int)blastresults[0]; //iinoka?
-
-// exit( 1 );
- }
- }
- else
- gappick0( mseq1[0], seq[scores->numinseq] );
- }
- else
- {
- table1 = (short *)calloc( tsize, sizeof( short ) );
- if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
- makecompositiontable_p( table1, scores[0].pointt );
- }
-
- selfscore0 = scores[0].selfscore;
- for( i=0; i<nin; i++ )
- {
- if( scores->orilen > scores[i].orilen )
- {
- longer = (double)scores->orilen;
- shorter = (double)scores[i].orilen;
- }
- else
- {
- longer = (double)scores[i].orilen; // nai
- shorter = (double)scores->orilen; //nai
- }
-
-#if LENFAC
- lenfac = 1.0 / ( shorter / longer * lenfacd + lenfacb / ( longer + lenfacc ) + lenfaca );
-// lenfac = 1.0 / ( (double)LENFACA + (double)LENFACB / ( (double)longer + (double)LENFACC ) + (double)shorter / (double)longer * LENFACD );
-// fprintf( stderr, "lenfac = %f l=%d,%d\n", lenfac,scores->orilen, scores[i].orilen );
-#else
- lenfac = 1.0;
-#endif
-
- if( doalign )
- {
- if( doalign == 'f' )
- {
- scores[i].score = ( 1.0 - blastresults[i] / MIN( scores->selfscore, scores[i].selfscore ) ) * 1;
- if( scores[i].score < 0.0 ) scores[i].score = 0.0;
- }
- else
- {
- if( fromaln )
- {
-// scores[i].score = ( 1.0 - (double)G__align11_noalign( amino_disLN, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[i].selfscore ) ) * 1;
- scores[i].score = ( 1.0 - (double)naivepairscore11( orialn[scores[i].numinseq], orialn[scores->numinseq], penalty ) / MIN( selfscore0, scores[i].selfscore ) ) * 1;
- }
- else
- {
- if( *depthpt == 0 ) fprintf( stderr, "\r%d / %d ", i, nin );
- gappick0( mseq2[0], seq[scores[i].numinseq] );
-// fprintf( stdout, "### before calc scores[%d] = %f (%c)\n", i, scores[i].score, qinoya == scores->numinseq?'o':'x' );
- scores[i].score = ( 1.0 - (double)G__align11_noalign( amino_disLN, -1200, -60, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[i].selfscore ) ) * 1;
-// fprintf( stderr, "scores[i] = %f\n", scores[i].score );
-// fprintf( stderr, "m1=%s\n", seq[scores[0].numinseq] );
-// fprintf( stderr, "m2=%s\n", seq[scores[i].numinseq] );
-// fprintf( stdout, "### before calc scores[%d] = %f (%c)\n", i, scores[i].score, qinoya == scores->numinseq?'o':'x' );
- }
- }
- }
- else
- {
- scores[i].score = ( 1.0 - (double)commonsextet_p( table1, scores[i].pointt ) / MIN( selfscore0, scores[i].selfscore ) ) * lenfac;
- if( scores[i].score > MAX6DIST ) scores[i].score = MAX6DIST;
- }
-// if( i ) fprintf( stderr, "%d-%d d %4.2f len %d %d\n", 1, i+1, scores[i].score, scores->orilen, scores[i].orilen );
- }
- if( doalign == 'f' ) free( blastresults );
- if( doalign == 0 ) free( table1 );
-//exit( 1 );
- }
-
-// fprintf( stderr, "sorting .. " );
- qsort( scores, nin, sizeof( Scores ), (int (*)())dcompare );
-// fprintf( stderr, "done.\n" );
-
-
- maxdist = scores[nin-1].score;
- if( fromaln ) // kanzen itch ga misalign sareteiru kamoshirenai.
- {
- if( scores[0].shimon == scores[nin-1].shimon && !strcmp( seq[scores[0].numinseq], seq[scores[nin-1].numinseq] ) )
- {
- maxdist = 0.0;
- }
-// fprintf( stderr, "maxdist?? = %f, nin=%d, %d inori\n", scores[nin-1].score, nin, scores[nin-1].numinseq+1 );
- }
-
-// fprintf( stderr, "maxdist? = %f, nin=%d\n", scores[nin-1].score, nin );
-
- if( nin == 1 ) fprintf( stderr, "nin=1, scores[0].score = %f\n", scores[0].score );
-
-// kokoni if( nin < 2 || ... )
-
- picks = AllocateIntVec( nin+1 );
- s_p_map = AllocateIntVec( nin+1 );
- s_y_map = AllocateIntVec( nin+1 );
- pickkouho = AllocateIntVec( nin+1 );
- closeh = AllocateIntVec( nin+1 );
-
-// nkouho = getkouho( pickkouho, (picksize+100)/nin, nin, scores, seq );
-// nkouho = getkouho( pickkouho, 1.0, nin, scores, seq ); // zenbu
-// fprintf( stderr, "selecting kouhos phase 2\n" );
-// if( nkouho == 0 )
-// {
-// fprintf( stderr, "selecting kouhos, phase 2\n" );
-// nkouho = getkouho( pickkouho, 1.0, nin, scores, seq );
-// }
-// fprintf( stderr, "\ndone\n\n" );
- for( i=0; i<nin; i++ ) pickkouho[i] = i+1; nkouho = nin-1; // zenbu
-
-
-
- iptr = picks;
- *iptr++ = 0;
- npick = 1;
- if( nkouho > 0 )
- {
-// fprintf( stderr, "pickkouho[0] = %d\n", pickkouho[0] );
-// fprintf( stderr, "pickkouho[nin-1] = %d\n", pickkouho[nin-1] );
- picktmp = pickkouho[nkouho-1];
-// fprintf( stderr, "\nMOST DISTANT kouho=%d, nin=%d, nkouho=%d\n", picktmp, nin, nkouho );
- nkouho--;
- if( ( scores[picktmp].shimon == scores[0].shimon ) && ( !strcmp( seq[scores[0].numinseq], seq[scores[picktmp].numinseq] ) ) )
- {
-// fprintf( stderr, "known, j=%d (%d inori)\n", 0, scores[picks[0]].numinseq );
-// fprintf( stderr, "%s\n%s\n", seq[scores[picktmp].numinseq], seq[scores[picks[0]].numinseq] );
- ;
- }
- else
- {
- *iptr++ = picktmp;
- npick++;
-// fprintf( stderr, "ok, %dth pick = %d (%d inori)\n", npick, picktmp, scores[picktmp].numinseq );
- }
- }
- i = 1;
- while( npick<picksize && nkouho>0 )
- {
- if( i )
- {
- i = 0;
- rn = nkouho * 0.5;
-// fprintf( stderr, "rn = %d\n", rn );
- }
- else
- {
- rn = rnd() * (nkouho);
- }
- picktmp = pickkouho[rn];
-// fprintf( stderr, "rn=%d/%d (%d inori), kouho=%d, nin=%d, nkouho=%d\n", rn, nkouho, scores[pickkouho[rn]].numinseq, pickkouho[rn], nin, nkouho );
-
-// fprintf( stderr, "#kouho before swap\n" );
-// for( i=0; i<nkouho; i++ ) fprintf( stderr, "%d ", pickkouho[i] ); fprintf( stderr, "\n" );
-
- nkouho--;
- pickkouho[rn] = pickkouho[nkouho];
-#if 1
-// fprintf( stderr, "#kouho after swap\n" );
-// for( i=0; i<nkouho; i++ ) fprintf( stderr, "%d ", pickkouho[i] ); fprintf( stderr, "\n" );
- for( j=0; j<npick; j++ )
- {
- if( scores[picktmp].shimon == scores[picks[j]].shimon && !strcmp( seq[scores[picks[j]].numinseq], seq[scores[picktmp].numinseq] ) )
- break;
- }
- if( j == npick )
-#endif
- {
-// fprintf( stderr, "ok, %dth pick = %d (%d inori)\n", npick, picktmp, scores[picktmp].numinseq );
- npick++;
- *iptr++ = picktmp;
- }
- else
- {
-// fprintf( stderr, "known, j=%d (%d inori)\n", j, scores[picks[j]].numinseq );
- }
- }
-#if 0
- for( i=0; i<nin; i++ )
- {
- fprintf( stderr, "i=%d/%d, scores[%d].score = %f, inori=%d\n", i, nin, i, scores[i].score, scores[i].numinseq );
- }
- fprintf( stderr, "range:nin=%d scores[%d].score <= %f\n", nin, npick, scores[nin-1].score);
- for( i=0; i<npick; i++ )
- {
- fprintf( stderr, "i=%d/%d, scores[%d].score = %f, inori=%d\n", i, npick, picks[i], scores[picks[i]].score, scores[picks[i]].numinseq );
- }
-exit( 1 );
-#endif
-
-// fprintf( stderr, "\nnkouho=%d, defaultq2 = %d (%d inori)\n", nkouho, picks[npick-1], scores[picks[npick-1]].numinseq );
-
- qsort( picks, npick, sizeof( int ), (int (*)())intcompare );
-
-// fprintf( stderr, "allocating..\n" );
-
-// fprintf( stderr, "allocating outs, npick = %d\n", npick );
- numin = calloc( npick, sizeof( int ) );
- tsukau = calloc( npick, sizeof( int ) );
- outs = calloc( npick, sizeof( Scores * ) );
- for( i=0; i<npick; i++ ) outs[i] = NULL;
- topol = AllocateIntCub( npick, 2, 0 );
- treeorder = AllocateIntVec( npick + 1 );
- len = AllocateFloatMtx( npick, 2 );
- pickmtx = AllocateFloatHalfMtx( npick );
-// yukomtx = AllocateFloatHalfMtx( npick );
-// minscoreinpick = AllocateDoubleVec( npick );
- yukos = AllocateIntVec( npick );
- p_o_map = AllocateIntVec( npick+1 );
- y_o_map = AllocateIntVec( npick+1 );
- hanni = AllocateFloatVec( npick );
-
- for( i=0; i<nin; i++ ) s_p_map[i] = -1;
-// fprintf( stderr, "npick = %d\n", npick );
-// fprintf( stderr, "picks =" );
- for( i=0; i<npick; i++ )
- {
- s_p_map[picks[i]] = i;
- p_o_map[i] = scores[picks[i]].numinseq;
-// fprintf( stderr, " %d (%dinori)\n", picks[i], scores[picks[i]].numinseq+1 );
- }
-// fprintf( stderr, "\n" );
-
-#if 0
- fprintf( stderr, "p_o_map =" );
- for( i=0; i<npick; i++ )
- {
- fprintf( stderr, " %d", p_o_map[i]+1 );
- }
- fprintf( stderr, "\n" );
- fprintf( stderr, "picks =" );
- for( i=0; i<npick; i++ )
- {
- fprintf( stderr, " %d", picks[i] );
- }
- fprintf( stderr, "\n" );
-#endif
-
- for( j=0; j<nin; j++ )
- {
- if( s_p_map[j] != -1 )
- {
- pickmtx[0][s_p_map[j]] = (float)scores[j].score;
-// fprintf( stderr, "pickmtx[0][%d] = %f\n", s_p_map[j], pickmtx[0][s_p_map[j]] );
- }
- }
-
- for( j=1; j<npick; j++ )
- {
- if( doalign )
- {
- if( doalign == 'f' )
- {
-// blastresults = callfasta( seq, scores, npick-j+1, picks+j-1, picks[j], 1 );
- blastresults = callfasta( seq, scores, npick, picks, picks[j], (j==1) );
- if( scores[picks[j]].selfscore != (int)blastresults[j] )
- {
- fprintf( stderr, "\n\nWARNING2: selfscore j=%d/%d\n", j, npick );
- fprintf( stderr, "scores[picks[j]].numinseq = %d\n", scores[picks[j]].numinseq+1 );
- fprintf( stderr, "scores[picks[j]].orilen = %d\n", scores[picks[j]].orilen );
- fprintf( stderr, "scores[picks[j]].selfscore = %d, but blastresults[j] = %f\n", scores[picks[j]].selfscore, blastresults[j] );
-// if( abs( scores[picks[j]].selfscore - (int)blastresults[j] ) > 2 )
-// exit( 1 );
-// scores->selfscore = (int)blastresults[0]; //iinoka?
- }
- }
- else
- gappick0( mseq1[0], seq[scores[picks[j]].numinseq] );
- }
- else
- {
- table1 = (short *)calloc( tsize, sizeof( short ) );
- if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
- makecompositiontable_p( table1, scores[picks[j]].pointt );
- }
-
- selfscore0 = scores[picks[j]].selfscore;
- pickmtx[j][0] = 0.0;
- for( i=j+1; i<npick; i++ )
- {
- if( scores[picks[j]].orilen > scores[picks[i]].orilen )
- {
- longer = (double)scores[picks[j]].orilen;
- shorter = (double)scores[picks[i]].orilen;
- }
- else
- {
- longer = (double)scores[picks[i]].orilen;
- shorter = (double)scores[picks[j]].orilen;
- }
-
- #if LENFAC
- lenfac = 1.0 / ( shorter / longer * lenfacd + lenfacb / ( longer + lenfacc ) + lenfaca );
- // lenfac = 1.0 / ( (double)LENFACA + (double)LENFACB / ( (double)longer + (double)LENFACC ) + (double)shorter / (double)longer * LENFACD );
- // fprintf( stderr, "lenfac = %f l=%d,%d\n", lenfac,scores->orilen, scores[i].orilen );
- #else
- lenfac = 1.0;
- #endif
-
- if( doalign )
- {
- if( doalign == 'f' )
- {
- pickmtx[j][i-j] = ( 1.0 - blastresults[i] / MIN( selfscore0, scores[picks[i]].selfscore ) ) * 1;
- if( pickmtx[j][i-j] < 0.0 ) pickmtx[j][i-j] = 0.0;
- }
- else
- {
- if( fromaln )
- {
- fprintf( stderr, "%d-%d/%d\r", j, i, npick );
- pickmtx[j][i-j] = ( 1.0 - (double)naivepairscore11( orialn[scores[picks[i]].numinseq], orialn[scores[picks[j]].numinseq], penalty ) / MIN( selfscore0, scores[picks[i]].selfscore ) ) * 1;
- }
- else
- {
-// fprintf( stderr, "\r%d / %d ", i, nin );
- gappick0( mseq2[0], seq[scores[picks[i]].numinseq] );
- pickmtx[j][i-j] = ( 1.0 - (double)G__align11_noalign( amino_disLN, -1200, -60, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[picks[i]].selfscore ) ) * 1;
- // fprintf( stderr, "scores[picks[i]] = %f\n", scores[picks[i]].score );
- }
- }
- }
- else
- {
- pickmtx[j][i-j] = ( 1.0 - (double)commonsextet_p( table1, scores[picks[i]].pointt ) / MIN( selfscore0, scores[picks[i]].selfscore ) ) * lenfac;
- if( pickmtx[j][i-j] > MAX6DIST ) pickmtx[j][i-j] = MAX6DIST;
- }
-
- }
- if( doalign == 'f' ) free( blastresults );
- if( doalign == 0 ) free( table1 );
- }
-
- dfromcp = AllocateDoubleMtx( npick, nin );
- dfromc = AllocateDoubleMtx( npick, 0 );
-
- for( i=0; i<npick; i++ ) for( j=0; j<nin; j++ )
- dfromcp[i][j] = -0.5;
- for( j=0; j<nin; j++ )
- {
- dfromcp[0][j] = ( scores[j].score );
-// fprintf( stderr, "j=%d, s_p_map[j]=%d\n", j, s_p_map[j] );
- }
-
- for( i=0; i<npick; i++ ) for( j=i; j<npick; j++ )
- {
- dfromcp[i][picks[j]] = dfromcp[j][picks[i]] = pickmtx[i][j-i];
- }
-
-#if 0
- fprintf( stderr, "pickmtx = \n" );
- for( i=0; i<npick; i++ )
- {
- for( j=i; j<npick; j++ )
- {
- fprintf( stderr, "pickmtx[%d][%d] = %f\n", p_o_map[i]+1, p_o_map[j]+1, pickmtx[i][j-i] );
- }
- }
- exit( 1 );
-#endif
-
-
-
-// for( i=0; i<npick-1; i++ ) for( j=i; j<npick; j++ )
-// fprintf( stderr, "dist[%d][%d] = %f\n", p_o_map[i]+1, p_o_map[j]+1, pickmtx[i][j-i] );
-
- for( i=0; i<npick; i++ ) tsukau[i] = 1;
- for( i=0; i<nin; i++ ) closeh[i] = -1;
- for( i=0; i<npick; i++ )
- {
- closeh[picks[i]] = picks[i];
-// fprintf( stderr, "i=%d/%d, picks[i]=%d, %d inori, closeh[%d] = %d \n", i, npick, picks[i], p_o_map[i]+1, picks[i], closeh[picks[i]] );
- }
-#if 0
- fprintf( stderr, "closeh = \n" );
- for( i=0; i<nin; i++ )
- {
- fprintf( stderr, "%d ", closeh[i] );
- }
- fprintf( stderr, "\n" );
-#endif
-#if DIANA
- for( i=0; i<npick-1; i++ ) for( j=i; j<npick; j++ )
- fprintf( stderr, "dist[%d][%d] = %f\n", p_o_map[i]+1, p_o_map[j]+1, pickmtx[i][j-i] );
- fprintf( stderr, "DIANA!!\n" );
- if( npick > 2 )
- {
- float avdist;
- float avdist1;
- float avdist2;
- float maxavdist;
- int splinter;
- int count;
- int dochokoho;
- splinter = 0;
- int *docholist;
- int *docholistbk;
- maxavdist = 0.0;
- for( i=0; i<npick; i++ )
- {
- avdist = 0.0;
- for( j=i+1; j<npick; j++ )
- {
- avdist += pickmtx[i][j-i];
- }
- for( j=0; j<i; j++ )
- {
- avdist += pickmtx[j][i-j];
- }
- avdist /= (npick-1);
- fprintf( stderr, "avdist[%d] = %f\n", p_o_map[i] + 1, avdist );
- if( maxavdist < avdist )
- {
- maxavdist = avdist;
- splinter = i;
- }
- }
- fprintf( stderr, "splinter = %d (%d inori), maxavdist = %f\n", splinter, p_o_map[splinter]+1, maxavdist );
-
- docholist = AllocateIntVec( npick );
- docholistbk = AllocateIntVec( npick );
- for( i=0; i<npick; i++ ) docholist[i] = 0;
- docholist[splinter] = 1;
- while( 1 )
- {
- for( i=0; i<npick; i++ ) docholistbk[i] = docholist[i];
- for( dochokoho = 0; dochokoho<npick; dochokoho++ )
- {
- fprintf( stderr, "dochokoho=%d\n", dochokoho );
- if( docholist[dochokoho] ) continue;
- count = 0;
- avdist1 = 0.0;
- i=dochokoho;
- {
- for( j=i+1; j<npick; j++ )
- {
- if( docholist[j] || j == dochokoho ) continue;
- avdist1 += pickmtx[i][j-i];
- count++;
- }
- for( j=0; j<i; j++ )
- {
- if( docholist[j] || j == dochokoho ) continue;
- avdist1 += pickmtx[j][i-j];
- count++;
- }
- }
- if( count < 1 ) avdist1 = 0.0;
- else avdist1 /= (float)count;
- fprintf( stderr, "docho %d (%dinori), avdist1 = %f\n", dochokoho, p_o_map[dochokoho] + 1, avdist1 );
-
- count = 0;
- avdist2 = 0.0;
- i=dochokoho;
- {
- for( j=i+1; j<npick; j++ )
- {
- if( !docholist[j] || j == dochokoho ) continue;
- avdist2 += pickmtx[i][j-i];
- count++;
- }
- for( j=0; j<i; j++ )
- {
- if( !docholist[j] || j == dochokoho ) continue;
- avdist2 += pickmtx[j][i-j];
- count++;
- }
- }
- if( count < 1 ) avdist2 = 0.0;
- else avdist2 /= (float)count;
- fprintf( stderr, "docho %d (%dinori), avdist2 = %f\n", dochokoho, p_o_map[dochokoho] + 1, avdist2 );
-
- if( avdist2 < avdist1 )
- {
- docholist[dochokoho] = 1;
- hanni[dochokoho] = avdist2;
- }
- else
- {
- docholist[dochokoho] = 0;
- hanni[dochokoho] = avdist1;
- }
- fprintf( stderr, "avdist1=%f, avdist2=%f\n", avdist1, avdist2 );
-
- }
- for( i=0; i<npick; i++ ) if( docholist[i] != docholistbk[i] ) break;
- if( i == npick ) break;
-
- fprintf( stderr, "docholist = \n" );
- for( i=0; i<npick; i++ ) fprintf( stderr, "%d ", docholist[i] );
- fprintf( stderr, "\n" );
- }
- fprintf( stderr, "docholist = \n" );
- for( i=0; i<npick; i++ ) fprintf( stderr, "%d ", docholist[i] );
- fprintf( stderr, "\n" );
-
- for( i=0; i<npick; i++ ) if( docholist[i] == 0 ) break;
- yukos[0] = picks[i];
- for( i=0; i<npick; i++ ) if( docholist[i] == 1 ) break;
- yukos[1] = picks[splinter];
-
- for( i=0; i<npick; i++ )
- {
- if( docholist[i] == 0 ) closeh[picks[i]] = yukos[0];
- if( docholist[i] == 1 ) closeh[picks[i]] = yukos[1];
- }
-// for( i=0; i<npick; i++ ) closeh[picks[i]] = -1; // CHUUI !! iminai
- nyuko = 2;
- free( docholist );
- free( docholistbk );
- }
- else if( npick > 1 )
- {
- nyuko = 2;
- yukos[0] = picks[0]; yukos[1] = picks[1];
- closeh[picks[0]] = yukos[0];
- closeh[picks[1]] = yukos[1];
- }
- else
- {
- nyuko = 1;
- yukos[0] = picks[0];
- closeh[picks[0]] = yukos[0];
- }
-#elif HUKINTOTREE
- if( npick > 2 )
- {
-#if 0
- float avdist;
- float maxavdist;
- int count;
- int splinter;
- maxavdist = 0.0;
- splinter=0;
- for( i=0; i<npick; i++ )
- {
- avdist = 0.0;
- for( j=i+1; j<npick; j++ )
- {
- avdist += pickmtx[i][j-i];
- }
- for( j=0; j<i; j++ )
- {
- avdist += pickmtx[j][i-j];
- }
- avdist /= (npick-1);
- fprintf( stderr, "avdist[%d] = %f\n", p_o_map[i] + 1, avdist );
- if( maxavdist < avdist )
- {
- maxavdist = avdist;
- splinter = i;
- }
- }
- fprintf( stderr, "splinter = %d (%d inori), maxavdist = %f\n", splinter, p_o_map[splinter]+1, maxavdist );
-#endif
-
-
-// fprintf( stderr, "check kaishi =>, npick=%d members = \n", npick );
-// for( i=0; i<npick; i++ ) fprintf( stderr, "%d (%d)", p_o_map[i]+1, picks[i] );
-// fprintf( stderr, "\n" );
- for( i=0; i<npick-1; i++ )
- {
- if( tsukau[i] == 0 ) continue;
- for( j=i+1; j<npick; j++ )
- {
-// float kijun = maxdist * 1/(npick-2);
-// float kijun = maxavdist * tokyoripara;
- float kijun;
- kijun = maxdist * tokyoripara; // atode kakunin
-// fprintf( stderr, "%d-%d\n", i, j );
-// fprintf( stderr, "maxdist = %f\n", maxdist );
-// if( i==0 && j == 1 ) continue; // machigai!! CHUUI!!
-// if( maxdist == pickmtx[i][j-i] ) continue;
- if( tsukau[j] == 0 ) continue;
-// fprintf( stderr, "checking %d-%d (%d-%d) %f, kijun=%f\n", p_o_map[i]+1, p_o_map[j]+1, i, j, pickmtx[i][j-i], kijun );
- if( pickmtx[i][j-i] < kijun )
- {
-// fprintf( stderr, "dame!! %d => %d, because %f < %f\n", p_o_map[j]+1, p_o_map[i]+1, pickmtx[i][j-i], kijun );
-#if 0
- if( scores[picks[i]].orilen > scores[picks[j]].orilen )
- {
- fprintf( stderr, "%d => %d\n", p_o_map[j]+1, p_o_map[i]+1 );
- tsukau[j] = 0;
- }
- else
- {
- fprintf( stderr, "%d => %d\n", p_o_map[i]+1, p_o_map[j]+1 );
- tsukau[i] = 0;
- }
- if( 0 && j == npick-1 ) tsukau[i] = 0;
- else tsukau[j] = 0;
- fprintf( stderr, "tsukau[%d] = %d (%d inori)\n", j, tsukau[j], p_o_map[j]+1 );
-#else
- tsukau[j] = 0;
- closeh[picks[j]] = closeh[picks[i]];
-// fprintf( stderr, "%d => tsukawanai\n", j );
-#endif
- }
- }
- }
- }
- for( ii=0,i=0; i<npick; i++ )
- {
- if( tsukau[i] )
- {
- dfromc[ii] = dfromcp[i];
- ii++;
- }
- else
- {
- free( dfromcp[i] );
- dfromcp[i] = NULL;
- }
- }
- dfromc[ii] = NULL;
-
- for( ii=0,i=0; i<npick; i++ )
- {
- if( tsukau[i] )
- {
- for( jj=ii,j=i; j<npick; j++ )
- {
- if( tsukau[j] )
- {
- pickmtx[ii][jj-ii] = pickmtx[i][j-i];
- jj++;
- }
- }
- ii++;
- }
- }
- for( ; ii<npick; ii++ )
- {
- free( pickmtx[ii] );
- pickmtx[ii] = NULL;
- }
-
- for( ii=0,i=0; i<npick; i++ )
- {
- if( tsukau[i] )
- {
- yukos[ii++] = picks[i];
- }
- }
-
-
- nyuko = ii;
- yukomtx = pickmtx;
- pickmtx = NULL;
-
-#endif
-#if 0
- for( i=0; i<npick; i++ ) for( j=i; j<npick; j++ )
- {
- if( tsukau[i] == 1 && tsukau[j] == 1 )
- fprintf( stderr, "dist[%d][%d] = %f (ok)\n", p_o_map[i]+1, p_o_map[j]+1, pickmtx[i][j-i] );
- else if( tsukau[i] == 0 && tsukau[j] == 0 )
- fprintf( stderr, "dist[%d][%d] = %f (xx)\n", p_o_map[i]+1, p_o_map[j]+1, pickmtx[i][j-i] );
- else
- fprintf( stderr, "%d-%d, okashii\n", p_o_map[i]+1, p_o_map[j]+1 );
- }
-#endif
-// FreeFloatHalfMtx( pickmtx, npick );
-
-
- for( i=0; i<nin; i++ ) s_y_map[i] = -1;
-// fprintf( stderr, "npick = %d\n", npick );
-// fprintf( stderr, "yukos =" );
- for( i=0; i<nyuko; i++ )
- {
- s_y_map[yukos[i]] = i;
- y_o_map[i] = scores[yukos[i]].numinseq;
-// fprintf( stderr, " %d\n", yukos[i] );
- }
-// fprintf( stderr, "\n" );
-#if 0
- for( i=0; i<nyuko; i++ )
- {
- fprintf( stderr, "y_o_map[%d] = %d\n", i, y_o_map[i]+1 );
- }
- for( i=0; i<nyuko; i++ ) for( j=i; j<nyuko; j++ )
- {
- fprintf( stderr, "yukodist[%d][%d] = %f (ok)\n", y_o_map[i]+1, y_o_map[j]+1, yukomtx[i][j-i] );
- }
-#endif
-
- for( i=0; i<nin; i++ )
- {
- if( closeh[i] != -1 )
- {
-// fprintf( stderr, "closeh[%d,%dinori] = %d,%dinori\n", i, scores[i].numinseq+1, closeh[i], scores[closeh[i]].numinseq+1 );
- }
- }
-
-#if 0
- for( i=0; i<nyuko; i++ )
- {
- minscoreinpick[i] = 99.9;
- for( j=i+1; j<nyuko; j++ )
- {
- if( minscoreinpick[i] > yukomtx[i][j-i] )
- minscoreinpick[i] = yukomtx[i][j-i];
- }
- for( j=0; j<i; j++ )
- {
- if( minscoreinpick[i] > yukomtx[j][i-j] )
- minscoreinpick[i] = yukomtx[j][i-j];
- }
- fprintf( stderr, "minscoreinpick[%d(%dinori)] = %f\n", i, y_o_map[i]+1, minscoreinpick[i] );
- }
-#endif
-
-
-#if TREE
- if( treeout )
- {
- children = calloc( nyuko+1, sizeof( char * ) );
- for( i=0; i<nyuko+1; i++ ) children[i] = NULL;
- }
-#endif
-// fprintf( stderr, "done..\n" );
-
-// fprintf( stderr, "classifying, nyuko=%d \n", nyuko );
- if( nyuko == 1 )
- {
- if( npick != 1 )
- {
- fprintf( stderr, "okashii, nyuko = 1, shikashi npick = %d\n", npick );
- exit( 1 );
- }
-// fprintf( stderr, "### itchi suru hazu, nazenara scores[nin-1].score=%f, selfscores=%d,%d\n", scores[nin-1].score, scores[nin-1].selfscore, scores->selfscore );
-// fprintf( stderr, "seq[%d] = scores->seq = \n%s\n", scores->numinseq, seq[scores->numinseq] );
-
- uniform = -1;
- for( j=0; j<nin; j++ )
- {
- belongto = 0;
- outs[belongto] = realloc( outs[belongto], sizeof( Scores ) * ( numin[belongto] + 1 ) );
- outs[belongto][numin[belongto]] = scores[j];
- numin[belongto]++;
- }
- }
- else
- {
-
-#if 0
- fprintf( stderr, "yukos = \n" );
- for( i=0; i<nyuko; i++ ) fprintf( stderr, "%d ", y_o_map[i] + 1 );
- fprintf( stderr, "\n" );
-#endif
- fprintf( stderr, "\n\n%dx%d distance matrix\n", nyuko, nin );
-
- for( i=1; i<nyuko; i++ )
- {
- fprintf( stderr, "%d / %d \r", i, nyuko );
-
- if( doalign )
- {
- if( doalign == 'f' )
- {
- blastresults = callfasta( seq, scores, nin, NULL, yukos[i], (i==1) );
- }
- else
- gappick0( mseq1[0], seq[scores[yukos[i]].numinseq] );
- }
- else
- {
- table1 = (short *)calloc( tsize, sizeof( short ) );
- if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
- makecompositiontable_p( table1, scores[yukos[i]].pointt );
- }
-
- selfscore0 = scores[yukos[i]].selfscore;
- for( j=0; j<nin; j++ )
- {
- if( scores[yukos[i]].orilen > scores[j].orilen )
- {
- longer = scores[yukos[i]].orilen;
- shorter = scores[j].orilen;
- }
- else
- {
- shorter = scores[yukos[i]].orilen;
- longer = scores[j].orilen;
- }
-
-#if LENFAC
-// lenfac = 1.0 / ( (double)LENFACA + (double)LENFACB / ( (double)longer + (double)LENFACC ) + (double)shorter / (double)longer * LENFACD );
- lenfac = 1.0 / ( shorter / longer * lenfacd + lenfacb / ( longer + lenfacc ) + lenfaca );
-// lenfac = 1.0 / ( shorter / longer * LENFACD + LENFACB / ( longer + LENFACC ) + LENFACA );
-// fprintf( stderr, "lenfac = %f, l=%d, %d\n", lenfac, scores[yukos[i]].orilen, scores[j].orilen );
-#else
- lenfac = 1.0;
-#endif
-#if 0 // iihazu -> dame
- ii = s_y_map[j]; jj=s_y_map[yukos[i]];
- if( ii != -1 && jj != -1 )
- {
- if( dfromc[ii][yukos[jj]] != -0.5 )
- {
- dfromc[i][j] = dfromc[ii][yukos[jj]];
- }
- else
- {
- if( ii > jj )
- {
- kk = jj;
- jj = ii;
- ii = kk;
- }
- dfromc[ii][yukos[jj]] =
- dfromc[i][j] = yukomtx[ii][jj-ii];
- }
- }
- else
-#else
- if( dfromc[i][j] == -0.5 )
-#endif
- {
- if( doalign )
- {
- if( doalign == 'f' )
- {
- dfromc[i][j] =
- ( 1.0 - blastresults[j] / MIN( selfscore0, scores[j].selfscore ) ) * 1;
- if( dfromc[i][j] < 0.0 ) dfromc[i][j] = 0.0;
- }
- else
- {
- if( fromaln )
- {
- dfromc[i][j] = ( 1.0 - (double)naivepairscore11( orialn[scores[j].numinseq], orialn[scores[yukos[i]].numinseq], penalty ) / MIN( selfscore0, scores[j].selfscore ) ) * 1;
- }
- else
- {
- gappick0( mseq2[0], seq[scores[j].numinseq] );
- dfromc[i][j] = ( 1.0 - (double)G__align11_noalign( amino_disLN, -1200, -60, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[j].selfscore ) ) * 1;
- }
- }
- }
- else
- {
- dfromc[i][j] = ( 1.0 - (double)commonsextet_p( table1, scores[j].pointt ) / MIN( selfscore0, scores[j].selfscore ) ) * lenfac;
- if( dfromc[i][j] > MAX6DIST ) dfromc[i][j] = MAX6DIST;
- }
- }
-// fprintf( stderr, "i,j=%d,%d (%d,%d)/ %d,%d, dfromc[][]=%f \n", i, j, scores[yukos[i]].numinseq+1, scores[j].numinseq+1, nyuko, nin, dfromc[i][j] );
-
-// if( i == 1 )
-// fprintf( stdout, "&&& dfromc[%d][%d] (%d,%d) = %f\n", i, j, p_o_map[i], scores[j].numinseq, dfromc[i][j] );
- }
-// fprintf( stderr, "i=%d, freeing\n", i );
- if( !doalign ) free( table1 );
- if( doalign && doalign == 'f' ) free( blastresults );
- }
- fprintf( stderr, " \r" );
-
-
-
-
- for( i=0; i<nyuko; i++ ) numin[i] = 0;
-// fprintf( stderr, "### itchi shinai hazu, nazenara scores[nin-1].score=%f, selfscores=%d,%d, len=%d,%d, nin=%d\n", scores[nin-1].score, scores[nin-1].selfscore, scores->selfscore, scores->orilen, scores[nin-1].orilen, nin );
- for( j=0; j<nin; j++ )
- {
-#if 0
- belongto = s_y_map[j];
- {
- fprintf( stderr, "belongto = %d (%dinori)\n", belongto, y_o_map[belongto]+1 );
- }
- if( belongto == -1 && closeh[j] != -1 )
-#endif
-#if 0
- if( closeh[j] != -1 )
- {
- belongto = s_y_map[closeh[j]];
-// if( belongto != -1 )
-// fprintf( stderr, "known, %d(%dinori)->%d(%dinori)\n", j, scores[j].numinseq+1, belongto, y_o_map[belongto]+1 );
- }
- else
-// if( belongto == -1 )
-#else
- belongto = s_y_map[j];
- if( belongto == -1 )
-#endif
- {
- belongto = 0; // default ha horyu
- minscore = dfromc[0][j];
- for( i=0; i<nyuko; i++ )
- {
-// fprintf( stderr, "checking %d/%d,%d/%d (%d-%d inori) minscore=%f, dfromc[0][j]=%f, dfromc[i][j]=%f\n", i, nyuko, j, nin, y_o_map[i], scores[j].numinseq, minscore, dfromc[0][j], dfromc[i][j] );
- if( scores[j].shimon == scores[yukos[i]].shimon && !strcmp( seq[scores[j].numinseq], seq[y_o_map[i]] ) )
- {
-// fprintf( stderr, "yuko-%d (%d in ori) to score-%d (%d inori) ha kanzen itch\n", i, y_o_map[i], j, scores[j].numinseq );
- belongto = i;
- break;
- }
- if( dfromc[i][j] < minscore )
-// if( dfromc[i][j] < minscore && minscore-dfromc[i][j] > ( minscoreinpick[yukos[i]] + minscoreinpick[j] ) * 1.0 )
-// if( rnd() < 0.5 ) // CHUUI !!!!!
- {
-// fprintf( stderr, "yuko-%d (%d in ori) to score-%d (%d inori) ha tikai, %f>%f\n", i, y_o_map[i]+1, j, scores[j].numinseq+1, minscore, dfromc[i][j] );
- minscore = dfromc[i][j];
- belongto = i;
- }
- }
- }
-#if 0
- if( dfromc[belongto][j] > minscoreinpick[belongto] )
- {
- fprintf( stderr, "dame, %f > %f\n", dfromc[belongto][j], minscoreinpick[belongto] );
- belongto = npick;
- }
- else
- fprintf( stderr, "ok, %f < %f\n", dfromc[belongto][j], minscoreinpick[belongto] );
-#endif
-// fprintf( stderr, "j=%d (%d inori) -> %d (%d inori) d=%f\n", j, scores[j].numinseq+1, belongto, y_o_map[belongto]+1, dfromc[belongto][j] );
-// fprintf( stderr, "numin = %d\n", numin[belongto] );
- outs[belongto] = realloc( outs[belongto], sizeof( Scores ) * ( numin[belongto] + 1 ) );
- outs[belongto][numin[belongto]] = scores[j];
- numin[belongto]++;
-
- }
- free( dfromcp );
- FreeDoubleMtx( dfromc );
-
-// fprintf( stderr, "##### npick = %d\n", npick );
-// fprintf( stderr, "##### nyuko = %d\n", nyuko );
-
-
- if( nyuko > 2 )
- {
- fprintf( stderr, "upgma " );
-// veryfastsupg_float_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len );
- fixed_musclesupg_float_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len, NULL );
- fprintf( stderr, "\r \r" );
- }
- else
- {
- topol[0][0] = (int *)realloc( topol[0][0], 2 * sizeof( int ) );
- topol[0][1] = (int *)realloc( topol[0][1], 2 * sizeof( int ) );
- topol[0][0][0] = 0;
- topol[0][0][1] = -1;
- topol[0][1][0] = 1;
- topol[0][1][1] = -1;
- }
- FreeFloatHalfMtx( yukomtx, npick );
-
-#if 0
- ii = nyuko-1;
- fprintf( stderr, "nyuko = %d, topol[][] = \n", nyuko );
- for( j=0; j<nyuko-1; j++ )
- {
- fprintf( stderr, "STEP%d \n", j );
- for( i=0; ; i++ )
- {
- fprintf( stderr, "%d ", ( topol[j][0][i] )+0 );
- if( topol[j][0][i] == -1 ) break;
- }
- fprintf( stderr, "\n" );
- for( i=0; ; i++ )
- {
- fprintf( stderr, "%d ", ( topol[j][1][i] )+0 );
- if( topol[j][1][i] == -1 ) break;
- }
- fprintf( stderr, "\n" );
- fprintf( stderr, "\n" );
- }
- exit( 1 );
-#endif
- jptr = treeorder;
- iptr = topol[nyuko-2][0]; while( *iptr != -1 ) *jptr++ = *iptr++;
- iptr = topol[nyuko-2][1]; while( *iptr != -1 ) *jptr++ = *iptr++;
- *jptr++ = -1;
- for( j=0; j<nyuko; j++ )
- {
-// fprintf( stderr, "treeorder[%d] = %d\n", j, treeorder[j] );
- if( treeorder[j] == -1 ) break;
- }
- }
-
- aligned = 1;
- for( i=0; i<nyuko; i++ )
- {
- ii = treeorder[i];
-#if 0
- if( numin[ii] > 1 )
- {
- fprintf( stderr, "\ncalling a child, pick%d (%d inori): # of mem=%d\n", i, p_o_map[ii]+1, numin[ii] );
- for( j=0; j<numin[ii]; j++ )
- {
- fprintf( stderr, "%d ", outs[ii][j].numinseq+1 );
- }
- fprintf( stderr, "\n" );
- }
-#endif
- aligned *= splitseq_mq( outs[ii], numin[ii], nlen, seq, orialn, name, inputfile, uniform, children+ii, alloclen, order, whichgroup, weight, depthpt, scores->numinseq );
- }
-
-
- for( i=0; i<nyuko; i++ )
- {
- if( !numin[i] )
- {
- fprintf( stderr, "i=%d/%d, ERROR!\n", i, nyuko );
- for( j=0; j<nyuko; j++ )
- fprintf( stderr, "numin[%d] = %d (rep=%d inori)\n", j, numin[j], y_o_map[j] );
- exit( 1 );
- }
- }
-
-#if TREE
- if( treeout )
- {
- treelen = 0;
- for( i=0; i<nyuko; i++ )
- treelen += strlen( children[i] );
- *tree = calloc( treelen + nin * 3, sizeof ( char ) );
- }
-#endif
-
-
- if( nin >= classsize || !aligned )
- val = 0;
- else
- val = 1;
-
- if( nyuko > 1 )
- {
- int *mem1p, *mem2p;
- int mem1size, mem2size;
- int v1 = 0, v2 = 0, v3 = 0;
- int nlim;
- int l;
- static int *mem1 = NULL;
- static int *mem2 = NULL;
- char **parttree = NULL; // by Mathog
-
-#if TREE
- if( treeout )
- {
- parttree = (char **)calloc( nyuko, sizeof( char * ) );
- for( i=0; i<nyuko; i++ )
- {
-// fprintf( stderr, "allocating parttree, size = %d\n", treelen + nin * 5 );
- parttree[i] = calloc( treelen + nin * 5, sizeof ( char ) );
- strcpy( parttree[i], children[i] );
- free( children[i] );
- }
- free( children );
- }
-#endif
- if( mem1 == NULL )
- {
- mem1 = AllocateIntVec( njob+1 );
- mem2 = AllocateIntVec( njob+1 );
- }
-
-// veryfastsupg_float_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len );
-
-// counteff_simple_float( nyuko, topol, len, eff );
-
-
- nlim = nyuko-1;
- for( l=0; l<nlim; l++ )
- {
- mem1p = topol[l][0];
- mptr = mem1;
- mem1size = 0;
- while( *mem1p != -1 )
- {
-// fprintf( stderr, "*mem1p = %d (%d inori), numin[]=%d\n", *mem1p, p_o_map[*mem1p], numin[*mem1p] );
- i = numin[*mem1p]; ptr = outs[*(mem1p++)];
- mem1size += i;
- while( i-- )
- {
- *mptr++ = (ptr++)->numinseq;
- }
- }
- *mptr = -1;
-
- mem2p = topol[l][1];
- mptr = mem2;
- mem2size = 0;
- while( *mem2p != -1 )
- {
-// fprintf( stderr, "*mem2p = %d (%d inori), numin[]=%d\n", *mem2p, p_o_map[*mem2p], numin[*mem2p] );
- i = numin[*mem2p]; ptr = outs[*(mem2p++)];
- mem2size += i;
- while( i-- )
- {
- *mptr++ = (ptr++)->numinseq;
- }
- }
- *mptr = -1;
-
- qsort( mem1, mem1size, sizeof( int ), (int (*)())intcompare );
- qsort( mem2, mem2size, sizeof( int ), (int (*)())intcompare );
-// selhead( mem1, numin[0] );
-// selhead( mem2, numin[1] );
-
-
-#if 0
- fprintf( stderr, "\n" );
- fprintf( stderr, "mem1 (nin=%d) = \n", nin );
- for( i=0; ; i++ )
- {
- fprintf( stderr, "%d ", mem1[i]+1 );
- if( mem1[i] == -1 ) break;
- }
- fprintf( stderr, "\n" );
- fprintf( stderr, "mem2 (nin=%d) = \n", nin );
- for( i=0; ; i++ )
- {
- fprintf( stderr, "%d ", mem2[i]+1 );
- if( mem2[i] == -1 ) break;
- }
- fprintf( stderr, "\n" );
-#endif
-
-#if 0
- fprintf( stderr, "before pairalign, l = %d, nyuko=%d, mem1size=%d, mem2size=%d\n", l, nyuko, mem1size, mem2size );
- fprintf( stderr, "before alignment\n" );
- for( j=0; j<mem1size; j++ )
- fprintf( stderr, "%s\n", seq[mem1[j]] );
- fprintf( stderr, "----\n" );
- for( j=0; j<mem2size; j++ )
- fprintf( stderr, "%s\n", seq[mem2[j]] );
- fprintf( stderr, "----\n\n" );
-#endif
-
- if( val )
- {
- fprintf( stderr, "\r Alignment %d-%d \r", mem1size, mem2size );
- if( *mem1 < *mem2 )
- pairalign( njob, nlen, seq, mem1, mem2, weight, alloclen );
- else
- pairalign( njob, nlen, seq, mem2, mem1, weight, alloclen );
- }
-
-#if TREE
- if( treeout )
- {
- v1 = topol[l][0][0];
- v2 = topol[l][1][0];
-
-// fprintf( stderr, "nyuko=%d, v1=%d, v2=%d\n", nyuko, v1, v2 );
- if( v1 > v2 )
- {
- v3 = v1;
- v1 = v2;
- v2 = v3;
- }
-// fprintf( stderr, "nyuko=%d, v1=%d, v2=%d after sort\n", nyuko, v1, v2 );
-// fprintf( stderr, "nyuko=%d, v1=%d, v2=%d\n", nyuko, v1, v2 );
-// fprintf( stderr, "v1=%d, v2=%d, parttree[v1]=%s, parttree[v2]=%s\n", v1, v2, parttree[v1], parttree[v2] );
- sprintf( *tree, "(%s,%s)", parttree[v1], parttree[v2] );
- strcpy( parttree[v1], *tree );
-// fprintf( stderr, "parttree[%d] = %s\n", v1, parttree[v1] );
-// fprintf( stderr, "*tree = %s\n", *tree );
- free( parttree[v2] ); parttree[v2] = NULL;
- }
-#endif
-
-#if 0
- fprintf( stderr, "after alignment\n" );
- for( j=0; j<mem1size; j++ )
- fprintf( stderr, "%s\n", seq[mem1[j]] );
- fprintf( stderr, "----\n" );
- for( j=0; j<mem2size; j++ )
- fprintf( stderr, "%s\n", seq[mem2[j]] );
- fprintf( stderr, "----\n\n" );
-#endif
- }
-#if TREE
- if( treeout )
- {
- free( parttree[v1] ); parttree[v1] = NULL;
-// fprintf( stderr, "*tree = %s\n", *tree );
-// FreeCharMtx( parttree );
- free( parttree ); parttree = NULL;
- }
-#endif
-
-#if 0
- fprintf( stderr, "after alignment\n" );
- for( i=0; i<nyuko; i++ )
- {
- for( j=0; j<numin[i]; j++ )
- fprintf( stderr, "%s\n", seq[outs[i][j].numinseq] );
- }
-#endif
-
- if( val )
- {
- groupid++;
- mptr = mem1; while( *mptr != -1 )
- {
-#if 0
- fprintf( stdout, "==g1-%d \n", *mptr+1 );
- fprintf( stdout, "%s \n", seq[*mptr] );
-#endif
- whichgroup[*mptr] = groupid;
- weight[*mptr++] *= 0.5;
- }
-
- mptr = mem2; while( *mptr != -1 )
- {
-#if 0
- fprintf( stdout, "=g2-%d ", *mptr+1 );
- fprintf( stdout, "%s \n", seq[*mptr] );
-#endif
- whichgroup[*mptr] = groupid;
- weight[*mptr++] *= 0.5;
- }
-
- if( numin[1] == 0 )
- {
- mptr = mem1; while( *mptr != -1 )
- {
- whichgroup[*mptr] = groupid;
- weight[*mptr++] *= (double)2.0/numin[0];
- }
- }
- }
- {
- if( *depthpt > maxdepth ) maxdepth = *depthpt;
- (*depthpt)++;
- }
- }
- else
- {
-#if TREE
- if( treeout )
- {
- sprintf( *tree, "%s", children[0] );
- free( children[0] );
- free( children );
- }
-#endif
- }
- for( i=0; i<npick; i++ ) free( (void *)outs[i] );
-// FreeFloatHalfMtx( pickmtx, npick );
-// FreeFloatHalfMtx( yukomtx, npick );
- FreeFloatMtx( len );
- FreeIntCub( topol );
- FreeIntVec( treeorder );
- free( outs );
- free( numin );
- free( picks );
- free( yukos );
- free( s_p_map );
- free( s_y_map );
- free( p_o_map );
- free( y_o_map );
- free( hanni );
- free( closeh );
- free( pickkouho );
- free( tsukau );
-// free( minscoreinpick );
- return val;
-}
-#endif
-
-static void alignparaphiles( int nseq, int *nlen, double *weight, char **seq, int nmem, int *members, int *alloclen )
-{
- int i, ilim;
- int *mem1 = AllocateIntVec( nmem );
- int *mem2 = AllocateIntVec( 2 );
-
- mem2[1] = -1;
- ilim = nmem-1;
- for( i=0; i<ilim; i++ )
- {
- mem1[i] = members[i];
- mem1[i+1] = -1;
- mem2[0] = members[i+1];
- pairalign( nseq, nlen, seq, mem1, mem2, weight, alloclen );
- }
- free( mem1 );
- free( mem2 );
-}
-
-
-
-
-
-
-
-
-
-
-
-int main( int argc, char *argv[] )
-{
- static char **name, **seq, **orialn;
- static int *grpseq;
- static char *tmpseq;
- static int **pointt;
- static int *nlen;
- int i, st, en;
- FILE *infp;
- FILE *treefp;
- char *treefile = NULL; //by Mathog
- char c;
- int alloclen;
- static int *order;
- static int *whichgroup;
- static double *weight;
- static char tmpname[B+100];
- int groupnum;
- int groupid;
- int pos;
- int *paramem;
- int npara;
- int completed;
- int orilen;
- int pscore;
- char *pt;
- int **tmpaminodis;
- static char com[1000];
- int depth;
- int aan;
-
- static Scores *scores;
- static short *table1;
- static char **tree;
-
-
-
- arguments( argc, argv );
-
- if( inputfile )
- {
- infp = fopen( inputfile, "r" );
- if( !infp )
- {
- fprintf( stderr, "Cannot open %s\n", inputfile );
- exit( 1 );
- }
- }
- else
- infp = stdin;
-
- getnumlen( infp );
- rewind( infp );
-
- if( njob < 2 )
- {
- fprintf( stderr, "At least 2 sequences should be input!\n"
- "Only %d sequence found.\n", njob );
- exit( 1 );
- }
-
-
- if( picksize == NOTSPECIFIED || picksize < 2 )
- picksize = PICKSIZE;
-
- if( classsize == NOTSPECIFIED || classsize < 0 )
- {
- classsize = njob + 1;
- }
- else
- {
-// picksize = MIN( picksize, (int)sqrt( classsize ) + 1);
- }
-
- if( tokyoripara == NOTSPECIFIED )
- {
- if( doalign ) tokyoripara = TOKYORIPARA_A;
- else tokyoripara = TOKYORIPARA;
- }
- if( picksize > njob )
- tokyoripara = 0.0;
-
-
- alloclen = nlenmax * 2;
- name = AllocateCharMtx( njob, B+1 );
-
- if( classsize == 1 )
- seq = AllocateCharMtx( njob, 0 );
- else
- seq = AllocateCharMtx( njob, alloclen+1 );
-
-
- nlen = AllocateIntVec( njob );
- tmpseq = calloc( nlenmax+1, sizeof( char ) );
- pointt = AllocateIntMtx( njob, 0 );
- grpseq = AllocateIntVec( nlenmax + 1 );
- order = (int *)calloc( njob + 1, sizeof( int ) );
- whichgroup = (int *)calloc( njob, sizeof( int ) );
- weight = (double *)calloc( njob, sizeof( double ) );
-
- fprintf( stderr, "alloclen = %d in main\n", alloclen );
-
- for( i=0; i<njob; i++ ) whichgroup[i] = 0;
- for( i=0; i<njob; i++ ) weight[i] = 1.0;
- for( i=0; i<njob; i++ ) order[i] = -1;
-
- if( classsize == 1 )
- readData_varlen( infp, name, nlen, seq );
- else
- readData_pointer( infp, name, nlen, seq );
-
- fclose( infp );
-
- if( fromaln ) doalign = 1;
-
- if( fromaln )
- {
- orialn = AllocateCharMtx( njob, alloclen+1 );
- for( i=0; i<njob; i++ )
- {
- if( strlen( seq[i] ) != nlenmax )
- {
- fprintf( stderr, "Input sequences must be aligned\n" );
- exit( 1 );
- }
- strcpy( orialn[i], seq[i] );
- }
- }
-
- constants( njob, seq );
-
- if( dorp == 'd' ) tsize = (int)pow( 4, 6 );
- else tsize = (int)pow( 6, 6 );
-
- if( dorp == 'd' )
- {
- lenfaca = DLENFACA;
- lenfacb = DLENFACB;
- lenfacc = DLENFACC;
- lenfacd = DLENFACD;
- }
- else
- {
- lenfaca = PLENFACA;
- lenfacb = PLENFACB;
- lenfacc = PLENFACC;
- lenfacd = PLENFACD;
- }
-
- maxl = 0;
- for( i=0; i<njob; i++ )
- {
- gappick0( tmpseq, seq[i] );
- nlen[i] = strlen( tmpseq );
- strcpy( seq[i], tmpseq );
- pointt[i] = AllocateIntVec( nlen[i]+1 ); // ??
-
- if( nlen[i] < 6 )
- {
- fprintf( stderr, "Seq %d, too short, %d characters\n", i+1, nlen[i] );
- fprintf( stderr, "name = %s\n", name[i] );
- fprintf( stderr, "seq = %s\n", seq[i] );
- exit( 1 );
-// continue;
- }
- if( nlen[i] > maxl ) maxl = nlen[i];
- if( dorp == 'd' ) /* nuc */
- {
- if( seq_grp_nuc( grpseq, tmpseq ) < 6 )
- {
- fprintf( stderr, "Seq %d, too short.\n", i+1 );
- fprintf( stderr, "name = %s\n", name[i] );
- fprintf( stderr, "seq = %s\n", seq[i] );
- exit( 1 );
-// continue;
- }
- makepointtable_nuc( pointt[i], grpseq );
- }
- else /* amino */
- {
- if( seq_grp( grpseq, tmpseq ) < 6 )
- {
- fprintf( stderr, "Seq %d, too short.\n", i+1 );
- fprintf( stderr, "name = %s\n", name[i] );
- fprintf( stderr, "seq = %s\n", seq[i] );
- exit( 1 );
-// continue;
- }
- makepointtable( pointt[i], grpseq );
- }
-// fprintf( stdout, ">%s\n", name[i] );
-// fprintf( stdout, "%s\n", seq[i] );
- }
-// exit( 1 );
-
-#if 0
- fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
-#endif
-
- initSignalSM();
-
- initFiles();
-
- WriteOptions( trap_g );
-
- c = seqcheck( seq );
- if( c )
- {
- fprintf( stderr, "Illeagal character %c\n", c );
- exit( 1 );
- }
-
- pid = (int)getpid();
- sprintf( datafile, "/tmp/data-%d", pid );
- sprintf( queryfile, "/tmp/query-%d", pid );
- sprintf( resultfile, "/tmp/fasta-%d", pid );
-
- scores = (Scores *)calloc( njob, sizeof( Scores ) );
-
-// fprintf( stderr, "\nCalculating i-i scores ... \n" );
- for( i=0; i<njob; i++ )
- {
- orilen = strlen( seq[i] );
- scores[i].numinseq = i; // irimasu
- scores[i].orilen = orilen;
- }
-
- if( doalign == 'f' )
- {
- fastapath = getenv( "FASTA_4_MAFFT" );
- if( !fastapath ) fastapath = "fasta34";
- fprintf( stderr, "fastapath=%s\n", fastapath );
- tmpaminodis = AllocateIntMtx( 0x80, 0x80 );
- getfastascoremtx( tmpaminodis );
- }
- else
- tmpaminodis = NULL;
-
- for( i=0; i<njob; i++ )
- {
- scores[i].pointt = pointt[i];
- scores[i].shimon = (int)pointt[i][0] + (int)pointt[i][1] + (int)pointt[i][scores[i].orilen-6];
- scores[i].name = name[i];
- if( doalign )
- {
- fprintf( stderr, "\r %05d/%05d ", i, njob );
- free( scores[i].pointt );
- if( doalign == 'f' )
- {
-#if 0
-#define KIZAMI 100
- int ipos = (int)( i / KIZAMI ) * KIZAMI;
- int iposamari = i % KIZAMI;
-
- fprintf( stderr, "%d / %d\r", i, njob );
-// fprintf( stderr, "ipos = %d\n", ipos );
-// fprintf( stderr, "iposamari = %d\n", iposamari );
-
-// fprintf( stderr, " calling blast, i=%d\n", i );
-// blastresults = callfasta( seq, scores+i, 1, 0, 1 );
- blastresults = callfasta( seq, scores+ipos, MIN(KIZAMI,njob-ipos), NULL, iposamari, (iposamari==0) );
-// fprintf( stderr, "done., i=%d\n\n", i );
- scores[i].selfscore = (int)blastresults[iposamari];
-#if 0
- for( j=0; j<100; j++ )
- {
- fprintf( stderr, "res[%d] = %f\n", j, blastresults[j] );
- }
-#endif
-// fprintf( stderr, "%d->selfscore = %d\n", i, scores[i].selfscore );
- free( blastresults );
-#else
- pscore = 0;
- if( scoremtx == -1 )
- {
- st = 1;
- en = 0;
- for( pt=seq[i]; *pt; pt++ )
- {
- if( *pt == 'u' ) *pt = 't';
- aan = amino_n[(int)*pt];
- if( aan<0 || aan >= 4 ) *pt = 'n';
-
- if( *pt == 'n' )
- {
- en++;
- if( st ) continue;
- else pscore += tmpaminodis[(int)*pt][(int)*pt];
- }
- else
- {
- st = 0;
- en = 0;
- pscore += tmpaminodis[(int)*pt][(int)*pt];
- }
- }
- scores[i].selfscore = pscore - en * tmpaminodis['n']['n'];
- }
- else
- {
- st = 1;
- en = 0;
- for( pt=seq[i]; *pt; pt++ )
- {
- aan = amino_n[(int)*pt];
- if( aan<0 || aan >= 20 ) *pt = 'X';
- if( *pt == 'X' )
- {
- en++;
- if( st ) continue;
- else pscore += tmpaminodis[(int)*pt][(int)*pt];
- }
- else
- {
- st = 0;
- en = 0;
- pscore += tmpaminodis[(int)*pt][(int)*pt];
- }
- }
- scores[i].selfscore = pscore - en * tmpaminodis['X']['X'];
- }
-#endif
- }
- else
- {
- pscore = 0;
- for( pt=seq[i]; *pt; pt++ )
- {
- pscore += amino_dis[(int)*pt][(int)*pt];
- }
- scores[i].selfscore = pscore;
- }
-// fprintf( stderr, "selfscore[%d] = %d\n", i+1, scores[i].selfscore );
- }
- else
- {
- table1 = (short *)calloc( tsize, sizeof( short ) );
- if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
- makecompositiontable_p( table1, pointt[i] );
- scores[i].selfscore = commonsextet_p( table1, pointt[i] );
- free( table1 );
- }
- }
- if( tmpaminodis ) FreeIntMtx( tmpaminodis );
-
- depth = 0;
-#if TREE
- if( treeout )
- {
- tree = (char **)calloc( 1, sizeof( char *) );
- *tree = NULL;
-// splitseq_bin( scores, njob, nlen, seq, name, inputfile, 0, tree, &alloclen, order, whichgroup, weight );
- completed = splitseq_mq( scores, njob, nlen, seq, orialn, name, inputfile, 0, tree, &alloclen, order, whichgroup, weight, &depth, -1 );
- treefile = (char *)calloc( strlen( inputfile ) + 10, sizeof( char ) );
- if( inputfile )
- sprintf( treefile, "%s.tree", inputfile );
- else
- sprintf( treefile, "splittbfast.tree" );
- treefp = fopen( treefile, "w" );
- fprintf( treefp, "%s\n", *tree );
- fclose( treefp );
- }
- else
- completed = splitseq_mq( scores, njob, nlen, seq, orialn, name, inputfile, 0, tree, &alloclen, order, whichgroup, weight, &depth, -1 );
-#else
- completed = splitseq_mq( scores, njob, nlen, seq, orialn, name, inputfile, 0, tree, &alloclen, order, whichgroup, weight, &depth, -1 );
-#endif
-
- fprintf( stderr, "\nDone.\n\n" );
-
-#if 1
- groupnum = 0;
- groupid = -1;
- paramem = NULL;
- npara = 0;
- for( i=0; i<njob; i++ )
- {
- pos = order[i];
- if( whichgroup[pos] != groupid )
- {
- groupnum++;
- groupid = whichgroup[pos];
- }
- if( whichgroup[pos] )
- {
- if( paramem )
- {
- paramem[npara] = -1;
- if( npara > 1 && classsize > 2 )
- {
- qsort( paramem, npara, sizeof( int ), (int (*)(const void *, const void*))intcompare );
-// selhead( paramem, npara );
- alignparaphiles( njob, nlen, weight, seq, npara, paramem, &alloclen );
- }
- free( paramem ); paramem = NULL; npara = 0;
- }
- sprintf( tmpname, "Group-%d %s", groupnum, name[pos]+1 );
- }
- else
- {
- paramem = realloc( paramem, sizeof( int) * ( npara + 2 ) );
- paramem[npara++] = pos;
- sprintf( tmpname, "Group-para %s", name[pos]+1 );
- }
- tmpname[B-1] = 0;
- if( classsize > 1 && classsize <= njob )
- strcpy( name[pos]+1, tmpname );
- }
- if( paramem )
- {
- paramem[npara] = -1;
- if( npara > 1 && classsize > 2 )
- {
- qsort( paramem, npara, sizeof( int ), (int (*)(const void *, const void*))intcompare );
-// selhead( paramem, npara );
- alignparaphiles( njob, nlen, weight, seq, npara, paramem, &alloclen );
- }
- free( paramem ); paramem = NULL; npara = 0;
- }
-#else
- for( i=0; i<njob; i++ )
- {
- sprintf( tmpname, "Group-%d %s", whichgroup[i], name[i]+1 );
- strcpy( name[i]+1, tmpname );
- }
-#endif
-
-
-// maketanni( name, seq, njob, nlenmax, nlen );
-
- fclose( trap_g );
-
-#if DEBUG
- fprintf( stderr, "writing alignment to stdout\n" );
-#endif
- if( reorder )
- writeData_reorder_pointer( stdout, njob, name, nlen, seq, order );
- else
- writeData_pointer( stdout, njob, name, nlen, seq );
-#if IODEBUG
- fprintf( stderr, "OSHIMAI\n" );
-#endif
- if( classsize == 1 )
- {
- fprintf( stderr, "\n\n" );
- fprintf( stderr, "----------------------------------------------------------------------------\n" );
- fprintf( stderr, "\n" );
- fprintf( stderr, "nseq = %d\n", njob );
- fprintf( stderr, "groupsize = %d, picksize=%d\n", classsize, picksize );
- fprintf( stderr, "The input sequences have been sorted so that similar sequences are close.\n" );
- if( reorder )
- fprintf( stderr, "The order of sequences has been changed according to estimated similarity.\n" );
-#if TREE
- if( treeout )
- {
- fprintf( stderr, "\n" );
- fprintf( stderr, "A guide tree is in the '%s' file.\n", treefile );
- }
-// else
-// {
-// fprintf( stderr, "To output guide tree,\n" );
-// fprintf( stderr, "%% %s -t -i %s\n", progName( argv[0] ), "inputfile" );
-// }
-#endif
- if( !doalign )
- {
- fprintf( stderr, "\n" );
- fprintf( stderr, "mafft --dpparttree might give a better result, although slow.\n" );
- fprintf( stderr, "mafft --fastaparttree is also available if you have fasta34.\n" );
- }
- fprintf( stderr, "\n" );
- fprintf( stderr, "----------------------------------------------------------------------------\n" );
- }
- else if( groupnum > 1 )
- {
- fprintf( stderr, "\n\n" );
- fprintf( stderr, "----------------------------------------------------------------------------\n" );
- fprintf( stderr, "\n" );
- fprintf( stderr, "groupsize = %d, picksize=%d\n", classsize, picksize );
- fprintf( stderr, "The input sequences have been classified into %d groups + some paraphyletic groups\n", groupnum );
- fprintf( stderr, "Note that the alignment is not completed.\n" );
- if( reorder )
- fprintf( stderr, "The order of sequences has been changed according to estimated similarity.\n" );
-#if TREE
- if( treeout )
- {
- fprintf( stderr, "\n" );
- fprintf( stderr, "A guide tree is in the '%s' file.\n", treefile );
- }
-// else
-// {
-// fprintf( stderr, "To output guide tree,\n" );
-// fprintf( stderr, "%% %s -t -i %s\n", progName( argv[0] ), "inputfile" );
-// }
-#endif
- if( !doalign )
- {
- fprintf( stderr, "\n" );
- fprintf( stderr, "mafft --dpparttree might give a better result, although slow.\n" );
- fprintf( stderr, "mafft --fastaparttree is also available if you have fasta34.\n" );
- }
- fprintf( stderr, "\n" );
- fprintf( stderr, "----------------------------------------------------------------------------\n" );
- }
- else
- {
- fprintf( stderr, "\n\n" );
- fprintf( stderr, "----------------------------------------------------------------------------\n" );
- fprintf( stderr, "\n" );
- fprintf( stderr, "nseq = %d\n", njob );
- fprintf( stderr, "groupsize = %d, partsize=%d\n", classsize, picksize );
-// fprintf( stderr, "A single alignment containing all the input sequences has been computed.\n" );
-// fprintf( stderr, "If the sequences are highly diverged and you feel there are too many gaps,\n" );
-// fprintf( stderr, "please try \n" );
-// fprintf( stderr, "%% mafft --parttree --groupsize 100 inputfile\n" );
-// fprintf( stderr, "which classifies the sequences into several groups with <~ 100 sequences\n" );
-// fprintf( stderr, "and performs only intra-group alignments.\n" );
- if( reorder )
- fprintf( stderr, "The order of sequences has been changed according to estimated similarity.\n" );
-#if TREE
- if( treeout )
- {
- fprintf( stderr, "\n" );
- fprintf( stderr, "A guide tree is in the '%s' file.\n", treefile );
- }
-// else
-// {
-// fprintf( stderr, "To output guide tree,\n" );
-// fprintf( stderr, "%% %s -t -i %s\n", progName( argv[0] ), "inputfile" );
-// }
-#endif
- if( !doalign || fromaln )
- {
- fprintf( stderr, "\n" );
- fprintf( stderr, "mafft --dpparttree might give a better result, although slow.\n" );
- fprintf( stderr, "mafft --fastaparttree is also available if you have fasta34.\n" );
- }
- fprintf( stderr, "\n" );
- fprintf( stderr, "----------------------------------------------------------------------------\n" );
- }
-#if TREE
- if( treeout ) free( treefile );
-#endif
-
-#if 0
- fprintf( stdout, "weight =\n" );
- for( i=0; i<njob; i++ )
- fprintf( stdout, "%d: %f\n", i+1, weight[i] );
-#endif
-
- if( doalign == 'f' )
- {
- strcpy( com, "rm -f" );
- strcat( com, " " );
- strcat( com, datafile );
- strcat( com, "* " );
- strcat( com, queryfile );
- strcat( com, " " );
- strcat( com, resultfile );
- fprintf( stderr, "%s\n", com );
- system( com );
- }
-
- SHOWVERSION;
-
- return( 0 );
-}
git://source.jalview.org / jabaws.git / blobdiff