+++ /dev/null
-#include "mltaln.h"
-
-#define TREE 1
-#define PICKSIZE 50 // must be >= 3
-#define WEIGHT 0
-#define HUKINTOTREE 1
-#define LENFAC 1
-
-// kouzoutai ni sasareru pointer ha static
-
-#define DEBUG 0
-#define IODEBUG 0
-#define SCOREOUT 0
-
-#define END_OF_VEC -1
-
-//static int doalign = 'f';
-static int doalign;
-static int treeout;
-static int classsize;
-static int picksize;
-static int maxl;
-static int tsize;
-static int numq;
-static int reorder;
-static int pid;
-static int maxdepth = 0;
-
-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 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;
- }
- }
-}
-
-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;
-}
-
-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 );
- 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 ) ErrorExit( "error in blast" );
-
- 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 );
-
-}
-
-static double *callfasta( 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 );
- }
-
-
- 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 );
-
- if( scoremtx == -1 )
- sprintf( com, "fasta34 -z3 -m10 -n -Q -b%d -E%d -d%d %s %s %d > %s\0", M, M, M, queryfile, datafile, 1, resultfile );
- else
- sprintf( com, "fastea34 -z3 -m10 -Q -b%d -E%d -d%d %s %s %d > %s\0", M, M, M, queryfile, datafile, 6, resultfile );
- res = system( com );
- if( res ) ErrorExit( "error in fasta" );
-
-//exit( 1 );
-
- 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 );
-}
-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 );
-}
-
-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;
-}
-
-void arguments( int argc, char *argv[] )
-{
- int c;
-
- doalign = 0;
- treeout = 0;
- 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;
- makedistmtx = 1;
- 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;
-
- 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;
-#if 0
- case 'm':
- fmodel = 1;
- break;
-#endif
- 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 'H':
- makedistmtx = 0;
- break;
- case 'e':
- fftscore = 0;
- break;
- case 'O':
- fftNoAnchStop = 1;
- break;
- case 'R':
- fftRepeatStop = 1;
- break;
- case 'Q':
- calledByXced = 1;
- break;
- case 'a':
- alg = 'a';
- break;
- case 'A':
- alg = 'A';
- break;
- case 'N':
- nevermemsave = 1;
- break;
- case 'M':
- alg = 'M';
- break;
- case 'B':
- break;
- case 'S':
- alg = 'S';
- 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 'z':
- fftThreshold = atoi( *++argv );
- --argc;
- goto nextoption;
- case 'w':
- fftWinSize = atoi( *++argv );
- --argc;
- goto nextoption;
- case 'Z':
- checkC = 1;
- break;
- 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;
-
-void seq_grp_nuc( int *grp, char *seq )
-{
- int tmp;
- while( *seq )
- {
- tmp = amino_grp[*seq++];
- if( tmp < 4 )
- *grp++ = tmp;
- else
- fprintf( stderr, "WARNING : Unknown character %c\r", *(seq-1) );
- }
- *grp = END_OF_VEC;
-}
-
-void seq_grp( int *grp, char *seq )
-{
- int tmp;
- while( *seq )
- {
- tmp = amino_grp[*seq++];
- if( tmp < 6 )
- *grp++ = tmp;
- else
- fprintf( stderr, "WARNING : Unknown character %c\r", *(seq-1) );
- }
- *grp = END_OF_VEC;
-}
-
-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 i;
- int l, len1, len2;
- int nlim;
- 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 ) FreeShortMtx( 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_noudp( 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 );
- }
- }
- 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 );
- break;
- case( 'A' ):
- if( clus1 == 1 && clus2 == 1 )
- {
-// fprintf( stderr, "%d-%d", clus1, clus2 );
- pscore = G__align11( mseq1, mseq2, *alloclen );
- }
- else
- {
-// fprintf( stderr, "%d-%d", clus1, clus2 );
- pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
- }
- 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
-
-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] );
- }
-}
-
-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 (3).\n" );
- else if( treemethod == 's' )
- fprintf( fp, "Tree = UPGMA (2).\n" );
- else if( treemethod == 'p' )
- fprintf( fp, "Tree = UPGMA (1).\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 **name, char *inputfile, int uniform, char **tree, int *alloclen, int *order, int *whichgroup, double *weight, int *depthpt )
-{
- int val;
- int ii, jj, kk;
- int *mptr;
- int treelen;
- int thisgroup;
- static int groupid = 0;
- static int branchid = 0;
- static int uniformid = 0;
- int i, j, k, jlim;
- int selfscore0, selfscore1;
- double **dfromc;
- float **pickmtx;
- float **pickmtx2;
- static short *table1;
- Scores **outs, *ptr;
- int *numin;
- int belongto;
- FILE *outfp;
- char *outputfile;
- char **children;
- char **parttree;
- char *tmptree;
- int *optr;
- static int *orderpos = NULL;
- int rn;
- int npick, npick2;
- int *picks;
- int *s_p_map;
- int *p_o_map;
- int nkouho;
- int *pickkouho;
- int *iptr;
- int *jptr;
- int aligned;
- int ***topol;
- int *treeorder;
- int picktmp;
- float **len;
- double minscore;
- double *minscoreinpick;
- double *hanni;
- int *saien;
- double *horyudist;
- double lenfac;
- double longer;
- double shorter;
- int off1, off2;
- char **mseq1, **mseq2;
- double *blastresults;
- static int palloclen = 0;
- double maxdist;
-
- if( orderpos == NULL )
- orderpos = order;
- if( palloclen == 0 )
- palloclen = *alloclen * 2;
-
-
-
- if( nin == 0 )
- {
-#if TREE
- if( treeout )
- {
- *tree = (char *)calloc( 1, sizeof( char ) );
- **tree = 0;
- }
-#endif
- return 1;
- }
- if( doalign && doalign != 'b' )
- {
- mseq1 = AllocateCharMtx( 1, palloclen );
- mseq2 = AllocateCharMtx( 1, palloclen );
- }
-
-
-
-#if 0
- okashii = 0;
- fprintf( stderr, "checking before swap!!\n" );
- for( i=0; i<nin; i++ )
- {
- fprintf( stderr, "scores[%d].seq (%d) = \n%s\n", i, scores[i].numinseq, scores[i].seq );
- if( strlen( seq[scores[i].numinseq] ) == 0 )
- {
- okashii = 1;
- fprintf( stderr, "OKASHII before swap!!\n" );
- }
- }
-#endif
-#if 0
- if( okashii == 1 )
- {
- for( i=0; i<njob; i++ )
- {
-// fprintf( stderr, "seq[%d] = \n%s\n", i, seq[i] );
- if( strlen( seq[i] ) == 0 )
- {
- fprintf( stderr, "OKASHII in seq!!\n" );
- }
- }
- exit( 1 );
- }
-#endif
-
- i = nin;
- ptr = scores;
- selfscore0 = scores->selfscore;
- belongto = 0;
- while( i-- )
- {
-// fprintf( stderr, "ptr-scores=%d, selfscore = %d, score = %f\n", ptr-scores, ptr->selfscore, 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
-
-
- if( doalign )
- {
- if( doalign == 'f' )
- {
- blastresults = callfasta( seq, scores, nin, 0, 1 );
- if( scores->selfscore != blastresults[0] )
- {
- fprintf( stderr, "WARNING: 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] );
-
- 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 == 'b' )
- {
- scores[i].score = ( 1.0 - blastresults[i] / MIN( scores->selfscore, scores[i].selfscore ) ) * 3;
- }
- else
- {
- fprintf( stderr, "\r%d / %d ", i, nin );
- gappick0( mseq2[0], seq[scores[i].numinseq] );
- scores[i].score = ( 1.0 - (double)L__align11_noalign( mseq1, mseq2, palloclen, &off1, &off2 ) / MIN( selfscore0, scores[i].selfscore ) ) * 3;
-// 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] );
- }
- }
- else
- scores[i].score = ( 1.0 - (double)commonsextet_p( table1, scores[i].pointt ) / MIN( selfscore0, scores[i].selfscore ) ) * lenfac;
-// 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 == 'b' ) 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;
- fprintf( stderr, "maxdist? = %f, nin=%d\n", scores[nin-1].score, nin );
-
- if( nin < 2 || uniform == -1 ) // kokodato chotto muda
- {
- 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 + 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;
- }
- picks = AllocateIntVec( nin+1 );
- s_p_map = AllocateIntVec( nin+1 );
- p_o_map = AllocateIntVec( nin+1 );
- pickkouho = 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
-
-
-// fprintf( stderr, "selecting picks" );
-
- 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] );
-// fprintf( stderr, "\nMOST DISTANT kouho=%d, nin=%d, nkouho=%d\n", pickkouho[nkouho], nin, nkouho );
- picktmp = pickkouho[nkouho-1];
- 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 );
-
- 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\n", picks[i] );
- }
-// fprintf( stderr, "\n" );
-
-#if 0
- fprintf( stderr, "p_o_map =" );
- for( i=0; i<npick; i++ )
- {
- fprintf( stderr, " %d", p_o_map[i] );
- }
- fprintf( stderr, "\n" );
-#endif
-// fprintf( stderr, "allocating..\n" );
-
-// fprintf( stderr, "allocating outs, npick = %d\n", npick );
- numin = calloc( npick+1, sizeof( int ) );
- outs = calloc( npick+1, sizeof( Scores * ) );
- for( i=0; i<npick+1; i++ ) outs[i] = NULL;
- topol = AllocateIntCub( npick+1, 2, 0 );
- treeorder = AllocateIntVec( npick+1 + 1 );
- len = AllocateFloatMtx( npick+1, 2 );
- pickmtx = AllocateFloatHalfMtx( npick+1 );
- minscoreinpick = AllocateDoubleVec( npick );
- hanni = AllocateDoubleVec( npick );
- saien = AllocateIntVec( npick );
- horyudist = AllocateDoubleVec( npick );
-#if TREE
- if( treeout )
- {
- children = calloc( npick+1, sizeof( char * ) );
- for( i=0; i<npick+1; i++ ) children[i] = NULL;
- }
-#endif
-// fprintf( stderr, "done..\n" );
-
-// fprintf( stderr, "classifying, pick=%d or %d \n", npick, npick+1 );
- if( npick == 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
- {
- dfromc = AllocateDoubleMtx( npick, nin );
-
- for( i=0; i<npick; i++ ) for( j=0; j<nin; j++ )
- dfromc[i][j] = -0.5;
- for( i=0; i<npick; i++ ) pickmtx[i][0] = 0.0;
- for( j=0; j<nin; j++ )
- {
- dfromc[0][j] = ( scores[j].score );
-// fprintf( stderr, "j=%d, s_p_map[j]=%d\n", j, s_p_map[j] );
- if( s_p_map[j] != -1 )
- {
- pickmtx[0][s_p_map[j]] = (float)dfromc[0][j];
- }
- }
-
- fprintf( stderr, "\n\n%dx%d distance matrix\n", npick, nin );
-
-#if 0
- fprintf( stderr, "picks = \n" );
- for( i=0; i<npick; i++ ) fprintf( stderr, "%d ", p_o_map[i] + 1 );
- fprintf( stderr, "\n" );
-
- fprintf( stderr, "scores = \n" );
- for( i=0; i<nin; i++ ) fprintf( stderr, "%d ", scores[i].numinseq + 1 );
- fprintf( stderr, "\n" );
-#endif
-
- for( i=1; i<npick; i++ )
- {
- fprintf( stderr, "%d / %d \r", i, npick );
-
- if( doalign )
- {
- if( doalign == 'f' )
- blastresults = callfasta( seq, scores, nin, picks[i], 0 );
- else
- gappick0( mseq1[0], seq[scores[picks[i]].numinseq] );
- }
- else
- {
- table1 = (short *)calloc( tsize, sizeof( short ) );
- if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
- makecompositiontable_p( table1, scores[picks[i]].pointt );
- }
-
- selfscore0 = scores[picks[i]].selfscore;
- for( j=0; j<nin; j++ )
- {
- if( scores[picks[i]].orilen > scores[j].orilen )
- {
- longer = scores[picks[i]].orilen;
- shorter = scores[j].orilen;
- }
- else
- {
- shorter = scores[picks[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[picks[i]].orilen, scores[j].orilen );
-#else
- lenfac = 1.0;
-#endif
-#if 1
- ii = s_p_map[j]; jj=s_p_map[picks[i]];
- if( ii != -1 && jj != -1 )
- {
- if( dfromc[ii][picks[jj]] != -0.5 )
- {
- dfromc[i][j] = dfromc[ii][picks[jj]];
- }
- else
- {
- if( doalign )
- {
- if( doalign == 'b' )
- {
- dfromc[ii][picks[jj]] =
- dfromc[i][j] =
- ( 1.0 - blastresults[j] / MIN( scores[picks[i]].selfscore, scores[j].selfscore ) ) * 3;
- }
- else
- {
- gappick0( mseq2[0], seq[scores[j].numinseq] );
- dfromc[ii][picks[jj]] =
- dfromc[i][j] =
- ( 1.0 - (double)L__align11_noalign( mseq1, mseq2, palloclen, &off1, &off2 ) / MIN( selfscore0, scores[j].selfscore ) ) * 3;
- }
- }
- else
- {
- dfromc[ii][picks[jj]] =
- dfromc[i][j] =
- ( 1.0 - (double)commonsextet_p( table1, scores[j].pointt ) / MIN( selfscore0, scores[j].selfscore ) ) * lenfac;
- }
- }
- }
- else
-#endif
- {
- if( doalign )
- {
- if( doalign == 'b' )
- {
- dfromc[i][j] =
- ( 1.0 - blastresults[j] / MIN( scores[picks[i]].selfscore, scores[j].selfscore ) ) * 3;
- }
- else
- {
- gappick0( mseq2[0], seq[scores[j].numinseq] );
- dfromc[i][j] = ( 1.0 - (double)L__align11_noalign( mseq1, mseq2, palloclen, &off1, &off2 ) / MIN( selfscore0, scores[j].selfscore ) ) * 3;
- }
- }
- else
- dfromc[i][j] = ( 1.0 - (double)commonsextet_p( table1, scores[j].pointt ) / MIN( selfscore0, scores[j].selfscore ) ) * lenfac;
- }
-// fprintf( stderr, "i,j=%d,%d (%d,%d)/ %d,%d, dfromc[][]=%f \n", i, j, scores[picks[i]].numinseq, scores[j].numinseq, npick, nin, dfromc[i][j] );
-
- if( s_p_map[j] != -1 )
- {
-// fprintf( stderr, "i=%d, j=%d, s_p_map[j]=%d, mtx[][]=%f\n", i, j, s_p_map[j], dfromc[i][j] );
- ii = i; jj=s_p_map[j];
- if( ii > jj )
- {
- kk = jj;
- jj = ii;
- ii = kk;
- }
- pickmtx[ii][jj-ii] = (float)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 == 'b' ) free( blastresults );
- }
- fprintf( stderr, " \r" );
-
- for( i=0; i<npick; i++ ) for( j=0; j<npick-i; j++ )
- fprintf( stderr, "pickmtx[%d][%d] = %f\n", p_o_map[i]+1, p_o_map[j+i]+1, pickmtx[i][j] );
-
- for( i=0; i<npick; i++ )
- {
- horyudist[i] = 0.0;
- hanni[i] = 0.00;
- saien[i] = -1;
- minscoreinpick[i] = 99.99;
- for( j=1; j<npick-i; j++ )
- {
-// fprintf( stderr, "minscoreinpick[%d] ?o %f\n", p_o_map[i]+1, pickmtx[i][j] );
- if( minscoreinpick[i] > pickmtx[i][j] )
- minscoreinpick[i] = pickmtx[i][j];
- }
- for( j=0; j<i; j++ )
- {
-// fprintf( stderr, "minscoreinpick[%d] ?x %f\n", p_o_map[i]+1, pickmtx[j][i-j] );
- if( minscoreinpick[i] > pickmtx[j][i-j] )
- minscoreinpick[i] = pickmtx[j][i-j];
- }
- minscoreinpick[i] *= 1.0;
- fprintf( stderr, "minscoreinpick[%d] = %f\n", p_o_map[i]+1, minscoreinpick[i] );
- }
-
-
-
-// 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++ )
- {
- belongto = s_p_map[j];
- if( belongto == -1 )
- {
- belongto = npick; // default ha horyu
- minscore = 9999.9;
- for( i=0; i<npick; i++ )
- {
-// fprintf( stderr, "checking %d/%d,%d/%d (%d-%d inori) minscore=%f, dfromc[0][j]=%f, dfromc[i][j]=%f\n", i, npick, j, nin, p_o_map[i], scores[j].numinseq, minscore, dfromc[0][j], dfromc[i][j] );
- if( scores[j].shimon == scores[picks[i]].shimon && !strcmp( seq[scores[j].numinseq], seq[p_o_map[i]] ) )
- {
-// fprintf( stderr, "pick-%d (%d in ori) to score-%d (%d inori) ha kanzen itch\n", i, p_o_map[i], j, scores[j].numinseq );
- belongto = i;
- break;
- }
- if( dfromc[i][j] < minscore )
-// if( rnd() < 0.5 ) // CHUUI !!!!!
- {
- if( dfromc[i][j] > hanni[i] )
- hanni[i] = dfromc[i][j];
- saien[i] = i;
- fprintf( stderr, "pick-%d (%d in ori) to score-%d (%d inori) ha tikai, %f>%f\n", i, p_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)\n", j, scores[j].numinseq+1, belongto, p_o_map[belongto]+1 );
- 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]++;
-
- }
-#if HUKINTOTREE
-
- for( i=0; i<npick; i++ )
- fprintf( stderr, "hanni[%d] = %f\n", p_o_map[i]+1, hanni[i] );
- 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++ )
- {
- fprintf( stderr, "numin[%d] = %d\n", p_o_map[i]+1, numin[i] );
- }
-
- if( npick > 2 )
- {
- npick2 = npick;
- for( i=0; i<npick-1; i++ )
- {
- if( numin[i] == 0 ) continue;
- for( j=i+1; j<npick; j++ )
- {
- fprintf( stderr, "%d-%d (%d-%d)\n", i, j, p_o_map[i]+1, p_o_map[j]+1 );
- if( i==0 && j == 1 )
- {
- fprintf( stderr, "continueing because i=0 and j=1\n" );
- continue;
- }
- if( numin[j] == 0 )
- {
- fprintf( stderr, "continueing because numin[j]=0\n" );
- continue;
- }
- fprintf( stderr, "checking %d-%d\n", p_o_map[i]+1, p_o_map[j]+1 );
- fprintf( stderr, "dist[i,j] = %f, hanni[i]=%f, hanni[j]=%f\n", pickmtx[i][j-i], hanni[i], hanni[j] );
-// if( pickmtx[i][j-i] < abs( dfromc[i][saien[i]] - dfromc[j][saien[i]] ) )
- if( pickmtx[i][j-i] < maxdist * 0.7 )
- {
- fprintf( stderr, "%d => %d, because %f<%f(md=%f)\n", p_o_map[j]+1, p_o_map[i]+1, pickmtx[i][j-i], maxdist*0.7, maxdist );
- npick2--;
- outs[i] = realloc( outs[i], sizeof( Scores ) * ( numin[i] + numin[j] ) );
- for( k=0; k<numin[j]; k++ ) outs[i][numin[i]+k] = outs[j][k];
- free( outs[j] ); outs[j] = NULL;
- numin[i] += numin[j];
- numin[j] = 0;
- }
-#if 0
- else if( pickmtx[i][j-i] < abs( dfromc[j][saien[j]] - dfromc[i][saien[i]] ) )
- {
- fprintf( stderr, "%d => %d\n", p_o_map[i]+1, p_o_map[j]+1 );
- npick2--;
- outs[j] = realloc( outs[j], sizeof( Scores ) * ( numin[j] + numin[i] ) );
- for( k=0; k<numin[i]; k++ ) outs[j][numin[j]+k] = outs[i][k];
- free( outs[i] ); outs[i] = NULL;
- numin[j] += numin[i];
- numin[i] = 0;
- }
-#endif
- else
- {
- ;
- }
- }
- }
- for( i=0; i<npick; i++ )
- {
- fprintf( stderr, "numin[%d] = %d, (%d inori)\n", i, numin[i], p_o_map[i]+1 );
- }
- pickmtx2 = AllocateFloatHalfMtx( npick2 );
- for( ii=0,i=0; i<npick; i++ )
- {
- if( numin[i] )
- {
- for( jj=ii,j=i; j<npick; j++ )
- {
- if( numin[j] )
- {
- pickmtx2[ii][jj-ii] = pickmtx[i][j-i];
- jj++;
- }
- }
- ii++;
- }
- }
- FreeFloatHalfMtx( pickmtx, npick );
- pickmtx = pickmtx2;
- for( ii=0,i=0; i<npick; i++ )
- {
- if( numin[i] )
- {
- outs[ii] = outs[i];
- numin[ii] = numin[i];
- p_o_map[ii] = p_o_map[i];
- ii++;
- }
- }
- for( ; ii<npick; ii++ )
- {
- outs[ii] = NULL;
- }
- for( i=0; i<npick; i++ )
- {
- fprintf( stderr, "numin[%d(%d)] = %d\n", i, p_o_map[i]+1, numin[i] );
- }
- npick = npick2;
-
- for( i=0; i<npick2; i++ ) for( j=i; j<npick2; j++ )
- fprintf( stderr, "dist[%d][%d] = %f\n", p_o_map[i]+1, p_o_map[j]+1, pickmtx[i][j-i] );
-
- }
-#endif
- FreeDoubleMtx( dfromc );
-
- fprintf( stderr, "##### npick = %d\n", npick );
-
-
- if( npick > 2 )
- {
- fprintf( stderr, "upgma " );
-// veryfastspg_float_realloc_nobk_halfmtx( npick, pickmtx, topol, len );
- veryfastsupg_float_realloc_nobk_halfmtx( npick, pickmtx, topol, len );
- 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;
-// free( pickmtx[0] );
- }
-
-#if 0
- ii = npick-1;
- fprintf( stderr, "npick = %d, topol[][] = \n", npick );
- for( j=0; j<npick; 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[npick-2][0]; while( *iptr != -1 ) *jptr++ = *iptr++;
- iptr = topol[npick-2][1]; while( *iptr != -1 ) *jptr++ = *iptr++;
- *jptr++ = -1;
- for( j=0; j<npick; j++ )
- {
-// fprintf( stderr, "treeorder[%d] = %d\n", j, treeorder[j] );
- if( treeorder[j] == -1 ) break;
- }
- }
-
- aligned = 1;
- for( i=0; i<npick; i++ )
- {
- ii = treeorder[i];
-#if 0
- fprintf( stderr, "\npick%d (%d inori): # of mem=%d\n", i, p_o_map[ii], numin[ii] );
- for( j=0; j<numin[ii]; j++ )
- {
- fprintf( stderr, "%d: %s\n", outs[ii][j].numinseq, seq[outs[ii][j].numinseq] );
- }
-#endif
- aligned *= splitseq_mq( outs[ii], numin[ii], nlen, seq, name, inputfile, uniform, children+ii, alloclen, order, whichgroup, weight, depthpt );
- }
-
-
- for( i=0; i<npick; i++ )
- {
- if( !numin[i] )
- {
- fprintf( stderr, "i=%d/%d, ERROR!\n", i, npick );
- for( j=0; j<npick; j++ )
- fprintf( stderr, "numin[%d] = %d (rep=%d inori)\n", j, numin[j], p_o_map[j] );
- exit( 1 );
- }
- }
-
-#if TREE
- if( treeout )
- {
- treelen = 0;
- for( i=0; i<npick; i++ )
- treelen += strlen( children[i] );
- *tree = calloc( treelen + nin * 3, sizeof ( char ) );
- }
-#endif
-
-
- if( nin >= classsize || !aligned )
- val = 0;
- else
- val = 1;
-
- if( npick > 1 )
- {
- int *mem1p, *mem2p;
- int mem1size, mem2size;
- int v1, v2, v3;
- int nlim;
- int l;
- int *mpt;
- static int *mem1 = NULL;
- static int *mem2 = NULL;
- char **parttree;
-
-#if TREE
- if( treeout )
- {
- parttree = (char **)calloc( npick, sizeof( char * ) );
- for( i=0; i<npick; 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( npick, pickmtx, topol, len );
-
-// counteff_simple_float( npick, topol, len, eff );
-
-
- nlim = npick-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, npick=%d, mem1size=%d, mem2size=%d\n", l, npick, 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
- fprintf( stderr, "\r Alignment %d-%d \r", mem1size, mem2size );
-
- if( val )
- {
- 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, "npick=%d, v1=%d, v2=%d\n", npick, v1, v2 );
- if( v1 > v2 )
- {
- v3 = v1;
- v1 = v2;
- v2 = v3;
- }
-// fprintf( stderr, "npick=%d, v1=%d, v2=%d after sort\n", npick, v1, v2 );
-// fprintf( stderr, "npick=%d, v1=%d, v2=%d\n", npick, 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<npick; 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 );
- FreeFloatMtx( len );
- FreeIntCub( topol );
- FreeIntVec( treeorder );
- free( outs );
- free( numin );
- free( picks );
- free( s_p_map );
- free( p_o_map );
- free( pickkouho );
- free( minscoreinpick );
- free( hanni );
- free( saien );
- free( horyudist );
- if( doalign && doalign != 'b' )
- {
- FreeCharMtx( mseq1 );
- FreeCharMtx( mseq2 );
- }
- return val;
-}
-#endif
-
-static void alignparaphiles( int nseq, int *nlen, double *weight, char **seq, int nmem, int *members, int *alloclen )
-{
- int i, j, 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;
- static int *grpseq;
- static char *tmpseq;
- static int **pointt;
- static int *nlen;
- int i, j;
- FILE *prep;
- FILE *infp;
- FILE *treefp;
- char *treefile;
- 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;
- double *blastresults;
- static char com[1000];
- int depth;
-
- 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 = PICKSIZE;
-
- if( classsize == NOTSPECIFIED || classsize == 0 )
- {
- classsize = njob + 1;
- }
- else
- {
-// picksize = MIN( picksize, (int)sqrt( classsize ) + 1);
- }
-
- alloclen = nlenmax * 2;
- name = AllocateCharMtx( njob, B+1 );
- seq = AllocateCharMtx( njob, alloclen+1 );
- nlen = AllocateIntVec( njob );
- tmpseq = calloc( nlenmax+1, sizeof( char ) );
- pointt = AllocateIntMtx( njob, nlenmax+1 );
- 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;
-
-// readData( infp, name, nlen, seq );
- readData_pointer( infp, name, nlen, seq );
-
- fclose( infp );
-
- 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 );
-
- 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 );
- }
- if( nlen[i] > maxl ) maxl = nlen[i];
- if( dorp == 'd' ) /* nuc */
- {
- seq_grp_nuc( grpseq, tmpseq );
- makepointtable_nuc( pointt[i], grpseq );
- }
- else /* amino */
- {
- seq_grp( grpseq, tmpseq );
- makepointtable( pointt[i], grpseq );
- }
- }
-
-#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\0", pid );
- sprintf( queryfile, "/tmp/query-%d\0", pid );
- sprintf( resultfile, "/tmp/fasta-%d\0", 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;
- }
-
- 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' )
- {
-#define KIZAMI 10
- 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), 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;
- 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 );
- }
- }
-
- 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, name, inputfile, 0, tree, &alloclen, order, whichgroup, weight, &depth );
- 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, name, inputfile, 0, tree, &alloclen, order, whichgroup, weight, &depth );
-#else
- completed = splitseq_mq( scores, njob, nlen, seq, name, inputfile, 0, tree, &alloclen, order, whichgroup, weight, &depth );
-#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;
- 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", njob );
- fprintf( stderr, "----------------------------------------------------------------------------\n", njob );
- fprintf( stderr, "\n", njob );
- 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 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
- fprintf( stderr, "\n", njob );
- fprintf( stderr, "----------------------------------------------------------------------------\n", njob );
- }
- else if( groupnum > 1 )
- {
- fprintf( stderr, "\n\n" );
- fprintf( stderr, "----------------------------------------------------------------------------\n", njob );
- 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 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
- fprintf( stderr, "\n" );
- fprintf( stderr, "----------------------------------------------------------------------------\n", njob );
- }
- else
- {
- fprintf( stderr, "\n\n" );
- fprintf( stderr, "----------------------------------------------------------------------------\n", njob );
- fprintf( stderr, "\n", njob );
- 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 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
- fprintf( stderr, "\n", njob );
- fprintf( stderr, "----------------------------------------------------------------------------\n", njob );
- }
-#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 == 'b' )
- {
- 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