+++ /dev/null
-#include "mltaln.h"
-#include "mtxutl.h"
-
-#define DEBUG 0
-#define TEST 0
-
-#define END_OF_VEC -1
-
-static int maxl;
-static int tsize;
-static char outputformat;
-static float 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
-
-void arguments( int argc, char *argv[] )
-{
- int c;
-
- inputfile = NULL;
- outputformat = 's';
- scoremtx = 1;
- nblosum = 62;
- dorp = NOTSPECIFIED;
- alg = 'X';
-
- while( --argc > 0 && (*++argv)[0] == '-' )
- {
- while ( (c = *++argv[0]) )
- {
- switch( c )
- {
- case 'i':
- inputfile = *++argv;
- fprintf( stderr, "inputfile = %s\n", inputfile );
- --argc;
- goto nextoption;
- case 'p':
- outputformat = 'p';
- break;
- case 'D':
- dorp = 'd';
- break;
- case 'P':
- dorp = 'p';
- break;
- default:
- fprintf( stderr, "illegal option %c\n", c );
- argc = 0;
- break;
- }
- }
- nextoption:
- ;
- }
- if( inputfile == NULL )
- {
- argc--;
- inputfile = *argv;
- fprintf( stderr, "inputfile = %s\n", inputfile );
- }
- if( argc != 0 )
- {
- fprintf( stderr, "Usage: mafft-distance [-PD] [-i inputfile] inputfile > outputfile\n" );
- exit( 1 );
- }
-}
-
-void 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;
- if( grp - grpbk < 6 )
- {
- fprintf( stderr, "\n\nWARNING: Too short.\nPlease also consider use mafft-ginsi, mafft-linsi or mafft-ginsi.\n\n\n" );
-// exit( 1 );
- *grpbk = -1;
- }
-}
-
-void 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;
- if( grp - grpbk < 6 )
- {
- fprintf( stderr, "\n\nWARNING: Too short.\nPlease also consider use mafft-ginsi, mafft-linsi or mafft-ginsi.\n\n\n" );
-// exit( 1 );
- *grpbk = -1;
- }
-}
-
-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( *pointt == -1 )
- return( 0 );
-
- 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;
-// fprintf( stderr, "cp - ct = %d (tsize = %d)\n", cp - ct, tsize );
- }
- *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;
-
- if( *n == -1 )
- {
- *pointt = -1;
- return;
- }
-
- 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;
-
- if( *n == -1 )
- {
- *pointt = -1;
- return;
- }
-
- 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;
-}
-
-int main( int argc, char **argv )
-{
- int i, j, initj;
- FILE *infp;
- char **seq;
- int *grpseq;
- char *tmpseq;
- int **pointt;
- static char **name;
- static int *nlen;
- double *mtxself;
- float score;
- static short *table1;
- float longer, shorter;
- float lenfac;
- float bunbo;
-
- arguments( argc, argv );
-
- if( inputfile )
- {
- infp = fopen( inputfile, "r" );
- if( !infp )
- {
- fprintf( stderr, "Cannot open %s\n", inputfile );
- exit( 1 );
- }
- }
- else
- infp = stdin;
-
-#if 0
- PreRead( stdin, &njob, &nlenmax );
-#else
- getnumlen( infp );
-#endif
- rewind( infp );
- if( njob < 2 )
- {
- fprintf( stderr, "At least 2 sequences should be input!\n"
- "Only %d sequence found.\n", njob );
- exit( 1 );
- }
-
- tmpseq = AllocateCharVec( nlenmax+1 );
- seq = AllocateCharMtx( njob, nlenmax+1 );
- grpseq = AllocateIntVec( nlenmax+1 );
- pointt = AllocateIntMtx( njob, nlenmax+1 );
- mtxself = AllocateDoubleVec( njob );
- pamN = NOTSPECIFIED;
- name = AllocateCharMtx( njob, B );
- nlen = AllocateIntVec( njob );
-
-#if 0
- FRead( infp, name, nlen, seq );
-#else
- readData_pointer( infp, name, nlen, seq );
-#endif
-
- 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 );
-// if( nlen[i] < 6 )
-// {
-// fprintf( stderr, "Seq %d, too short, %d characters\n", i+1, nlen[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 );
- }
- }
- fprintf( stderr, "\nCalculating i-i scores ... " );
- for( i=0; i<njob; i++ )
- {
- table1 = (short *)calloc( tsize, sizeof( short ) );
- if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
- makecompositiontable_p( table1, pointt[i] );
-
- score = commonsextet_p( table1, pointt[i] );
- mtxself[i] = score;
- free( table1 );
- }
-
- fprintf( stderr, "done.\n" );
- fprintf( stderr, "\nCalculating i-j scores ... \n" );
- if( outputformat == 'p' ) fprintf( stdout, "%-5d", njob );
- for( i=0; i<njob; i++ )
- {
- if( outputformat == 'p' ) fprintf( stdout, "\n%-9d ", i+1 );
- table1 = (short *)calloc( tsize, sizeof( short ) );
- if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
- if( i % 10 == 0 )
- {
- fprintf( stderr, "%4d / %4d\r", i+1, njob );
- }
- makecompositiontable_p( table1, pointt[i] );
-
-
- if( outputformat == 'p' ) initj = 0;
- else initj = i+1;
- for( j=initj; j<njob; j++ )
- {
- if( nlen[i] > nlen[j] )
- {
- longer=(float)nlen[i];
- shorter=(float)nlen[j];
- }
- else
- {
- longer=(float)nlen[j];
- shorter=(float)nlen[i];
- }
-// lenfac = 3.0 / ( LENFACA + LENFACB / ( longer + LENFACC ) + shorter / longer * LENFACD );
- lenfac = 1.0 / ( shorter / longer * lenfacd + lenfacb / ( longer + lenfacc ) + lenfaca );
-// lenfac = 1.0;
-// fprintf( stderr, "lenfac = %f (%.0f,%.0f)\n", lenfac, longer, shorter );
- score = commonsextet_p( table1, pointt[j] );
- bunbo = MIN( mtxself[i], mtxself[j] );
- if( outputformat == 'p' )
- {
- if( bunbo == 0.0 )
- fprintf( stdout, " %8.6f", 1.0 );
- else
- fprintf( stdout, " %8.6f", ( 1.0 - score / bunbo ) * lenfac );
- if( j % 7 == 6 ) fprintf( stdout, "\n" );
- }
- else
- {
- if( bunbo == 0.0 )
- fprintf( stdout, "%d-%d d=%4.2f l=%d,%d\n", i+1, j+1, 1.0, nlen[i], nlen[j] );
- else
- fprintf( stdout, "%d-%d d=%4.2f l=%d,%d\n", i+1, j+1, ( 1.0 - score / bunbo ) * lenfac, nlen[i], nlen[j] );
- }
-// fprintf( stderr, "##### mtx = %f, mtx[i][0]=%f, mtx[j][0]=%f, bunbo=%f\n", mtx[i][j-i], mtx[i][0], mtx[j][0], bunbo );
-// score = (double)commonsextet_p( table1, pointt[j] );
-// fprintf( stdout, "%d-%d d=%4.2f l=%d,%d\n", i+1, j+1, ( 1.0 - score / MIN( mtxself[i], mtxself[j] ) ) * 3, nlen[i], nlen[j] );
-
-
- }
- free( table1 );
- }
-
- fprintf( stderr, "\n" );
- if( outputformat == 'p' ) fprintf( stdout, "\n" );
- SHOWVERSION;
- exit( 0 );
-}
git://source.jalview.org / jabaws.git / blobdiff