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
}
}
-#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 )
{
tbrweight = 3;
checkC = 0;
treemethod = 'X';
+ sueff_global = 0.1;
contin = 0;
scoremtx = 1;
kobetsubunkatsu = 0;
dorp = NOTSPECIFIED;
ppenalty = -1530;
ppenalty_ex = NOTSPECIFIED;
+ penalty_shift_factor = 1000.0;
poffset = -123;
kimuraR = NOTSPECIFIED;
pamN = NOTSPECIFIED;
classsize = NOTSPECIFIED;
picksize = NOTSPECIFIED;
tokyoripara = NOTSPECIFIED;
+ legacygapcost = 0;
+ nwildcard = 0;
+ outnumber = 0;
while( --argc > 0 && (*++argv)[0] == '-' )
{
switch( c )
{
case 'p':
- picksize = atoi( *++argv );
+ picksize = myatoi( *++argv );
fprintf( stderr, "picksize = %d\n", picksize );
--argc;
goto nextoption;
case 's':
- classsize = atoi( *++argv );
+ classsize = myatoi( *++argv );
fprintf( stderr, "groupsize = %d\n", classsize );
--argc;
goto nextoption;
// fprintf( stderr, "ppenalty = %d\n", ppenalty );
--argc;
goto nextoption;
+ case 'Q':
+ penalty_shift_factor = atof( *++argv );
+ --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 'k':
- kimuraR = atoi( *++argv );
+ kimuraR = myatoi( *++argv );
fprintf( stderr, "kimuraR = %d\n", kimuraR );
--argc;
goto nextoption;
case 'b':
- nblosum = atoi( *++argv );
+ nblosum = myatoi( *++argv );
scoremtx = 1;
// fprintf( stderr, "blosum %d\n", nblosum );
--argc;
goto nextoption;
case 'j':
- pamN = atoi( *++argv );
+ pamN = myatoi( *++argv );
scoremtx = 0;
TMorJTT = JTT;
fprintf( stderr, "jtt %d\n", pamN );
--argc;
goto nextoption;
case 'm':
- pamN = atoi( *++argv );
+ pamN = myatoi( *++argv );
scoremtx = 0;
TMorJTT = TM;
fprintf( stderr, "tm %d\n", pamN );
case 'l':
uselongest = 0;
break;
+ case 'n' :
+ outnumber = 1;
+ break;
#if 1
case 'a':
fmodel = 1;
case 'Z':
fromaln = 1;
break;
- case 'L':
+ case 'U':
doalign = 1;
break;
case 'x':
case 'O':
fftNoAnchStop = 1;
break;
+ case 'L':
+ legacygapcost = 1;
+ break;
#if 0
case 'R':
fftRepeatStop = 1;
case 'a':
alg = 'a';
break;
-#endif
case 'R':
alg = 'R';
break;
case 'Q':
alg = 'Q';
break;
+#endif
case 'A':
alg = 'A';
break;
tbutree = 0;
break;
case 'X':
- treemethod = 'X'; // mix
- break;
+ treemethod = 'X'; // tsukawareteiru ????
+ sueff_global = atof( *++argv );
+ fprintf( stderr, "sueff_global = %f\n", sueff_global );
+ --argc;
+ goto nextoption;
case 'E':
treemethod = 'E'; // upg (average)
break;
treemethod = 'q'; // minimum
break;
case 'z':
- fftThreshold = atoi( *++argv );
+ fftThreshold = myatoi( *++argv );
--argc;
goto nextoption;
case 'w':
- fftWinSize = atoi( *++argv );
+ fftWinSize = myatoi( *++argv );
--argc;
goto nextoption;
+ case ':':
+ nwildcard = 1;
+ break;
default:
fprintf( stderr, "illegal option %c\n", c );
argc = 0;
}
}
-static int maxl;
-static int tsize;
+static int nunknown = 0;
int seq_grp_nuc( int *grp, char *seq )
{
if( tmp < 4 )
*grp++ = tmp;
else
- fprintf( stderr, "WARNING : Unknown character %c\r", *(seq-1) );
+ nunknown++;
}
*grp = END_OF_VEC;
return( grp-grpbk );
if( tmp < 6 )
*grp++ = tmp;
else
- fprintf( stderr, "WARNING : Unknown character %c\r", *(seq-1) );
+ nunknown++;
}
*grp = END_OF_VEC;
return( grp-grpbk );
table[point]++;
}
-int commonsextet_p( short *table, int *pointt )
+static int localcommonsextet_p( short *table, int *pointt )
{
int value = 0;
short tmp;
{
int l, len1, len2;
int clus1, clus2;
- float pscore, tscore;
+ double 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;
+// double dumfl = 0.0;
+ double dumdb = 0.0;
int ffttry;
int m1, m2;
#if 0
}
#if WEIGHT
- clus1 = fastconjuction_noname( mem1, seq, mseq1, effarr1, weight, indication1 );
- clus2 = fastconjuction_noname( mem2, seq, mseq2, effarr2, weight, indication2 );
+ clus1 = fastconjuction_noname( mem1, seq, mseq1, effarr1, weight, indication1, 0.0 );
+ clus2 = fastconjuction_noname( mem2, seq, mseq2, effarr2, weight, indication2, 0.0 );
#else
clus1 = fastconjuction_noweight( mem1, seq, mseq1, effarr1, indication1 );
clus2 = fastconjuction_noweight( mem2, seq, mseq2, effarr2, indication2 );
fprintf( stderr, "\nlen1=%d, len2=%d, Switching to the memsave mode\n", len1, len2 );
alg = 'M';
if( commonIP ) FreeIntMtx( commonIP );
+ commonIP = 0;
commonAlloc1 = 0;
commonAlloc2 = 0;
}
{
fprintf( stderr, "\bm" );
// fprintf( stderr, "%d-%d", clus1, clus2 );
- pscore = Falign_udpari_long( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1 );
+ pscore = Falign_udpari_long( NULL, NULL, n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, NULL, NULL, 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 );
+ pscore = Falign( NULL, NULL, n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, NULL, NULL, clus1, clus2, *alloclen, fftlog+m1, NULL, 0, NULL );
}
}
else
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 );
- }
+ pscore = MSalignmm( n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
break;
case( 'A' ):
if( clus1 == 1 && clus2 == 1 )
{
// fprintf( stderr, "%d-%d", clus1, clus2 );
- pscore = G__align11( mseq1, mseq2, *alloclen, outgap, outgap );
+ pscore = G__align11( n_dis_consweight_multi, 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 );
+ pscore = A__align( n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumdb, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap, -1, -1 );
}
break;
default:
int selfscore0;
double **dfromc;
double **dfromcp;
- float **pickmtx;
- float **yukomtx;
+ double **pickmtx;
+ double **yukomtx;
static short *table1;
Scores **outs, *ptr;
int *numin;
int ***topol;
int *treeorder;
int picktmp;
- float **len;
+ double **len;
double minscore;
// double *minscoreinpick;
- float *hanni;
+ double *hanni;
double lenfac;
double longer;
double shorter;
static char **mseq2 = NULL;
double *blastresults = NULL; // by Mathog, a guess
static int palloclen = 0;
- float maxdist;
+ double maxdist;
if( orderpos == NULL )
orderpos = order;
}
free( tmptree );
- *tree = (char *)calloc( treelen + nin + 5, sizeof( char ) );
- if( nin > 1 ) **tree = '(';
- else **tree = '\0';
-// **tree = '\0';
+ *tree = (char *)calloc( treelen + nin + 15, sizeof( char ) );
+ **tree = '\n';
+ if( nin > 1 )
+ {
+ *(*tree+1) = '(';
+ *(*tree+2) = '\0';
+ }
+ else
+ {
+ *(*tree+1) = '\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, ")" );
+ if( nin > 1 ) strcat( *tree, ")\n" );
+ else strcat( *tree, "\n" );
// fprintf( stdout, "*tree = %s\n", *tree );
}
{
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)G__align11_noalign( n_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;
+ scores[i].score = ( 1.0 - (double)G__align11_noalign( n_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] );
}
else
{
- scores[i].score = ( 1.0 - (double)commonsextet_p( table1, scores[i].pointt ) / MIN( selfscore0, scores[i].selfscore ) ) * lenfac;
+ scores[i].score = ( 1.0 - (double)localcommonsextet_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( s_p_map[j] != -1 )
{
- pickmtx[0][s_p_map[j]] = (float)scores[j].score;
+ pickmtx[0][s_p_map[j]] = (double)scores[j].score;
// fprintf( stderr, "pickmtx[0][%d] = %f\n", s_p_map[j], pickmtx[0][s_p_map[j]] );
}
}
{
// 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;
+ pickmtx[j][i-j] = ( 1.0 - (double)G__align11_noalign( n_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;
+ pickmtx[j][i-j] = ( 1.0 - (double)localcommonsextet_p( table1, scores[picks[i]].pointt ) / MIN( selfscore0, scores[picks[i]].selfscore ) ) * lenfac;
if( pickmtx[j][i-j] > MAX6DIST ) pickmtx[j][i-j] = MAX6DIST;
}
fprintf( stderr, "DIANA!!\n" );
if( npick > 2 )
{
- float avdist;
- float avdist1;
- float avdist2;
- float maxavdist;
+ double avdist;
+ double avdist1;
+ double avdist2;
+ double maxavdist;
int splinter;
int count;
int dochokoho;
}
}
if( count < 1 ) avdist1 = 0.0;
- else avdist1 /= (float)count;
+ else avdist1 /= (double)count;
fprintf( stderr, "docho %d (%dinori), avdist1 = %f\n", dochokoho, p_o_map[dochokoho] + 1, avdist1 );
count = 0;
}
}
if( count < 1 ) avdist2 = 0.0;
- else avdist2 /= (float)count;
+ else avdist2 /= (double)count;
fprintf( stderr, "docho %d (%dinori), avdist2 = %f\n", dochokoho, p_o_map[dochokoho] + 1, avdist2 );
if( avdist2 < avdist1 )
if( npick > 2 )
{
#if 0
- float avdist;
- float maxavdist;
+ double avdist;
+ double maxavdist;
int count;
int splinter;
maxavdist = 0.0;
if( tsukau[i] == 0 ) continue;
for( j=i+1; j<npick; j++ )
{
-// float kijun = maxdist * 1/(npick-2);
-// float kijun = maxavdist * tokyoripara;
- float kijun;
+// double kijun = maxdist * 1/(npick-2);
+// double kijun = maxavdist * tokyoripara;
+ double kijun;
kijun = maxdist * tokyoripara; // atode kakunin
// fprintf( stderr, "%d-%d\n", i, j );
// fprintf( stderr, "maxdist = %f\n", maxdist );
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;
+ dfromc[i][j] = ( 1.0 - (double)G__align11_noalign( n_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;
+ dfromc[i][j] = ( 1.0 - (double)localcommonsextet_p( table1, scores[j].pointt ) / MIN( selfscore0, scores[j].selfscore ) ) * lenfac;
if( dfromc[i][j] > MAX6DIST ) dfromc[i][j] = MAX6DIST;
}
}
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 );
+// veryfastsupg_double_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len );
+ fixed_musclesupg_double_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len, NULL, 1, 1 );
fprintf( stderr, "\r \r" );
}
else
mem2 = AllocateIntVec( njob+1 );
}
-// veryfastsupg_float_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len );
+// veryfastsupg_double_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len );
-// counteff_simple_float( nyuko, topol, len, eff );
+// counteff_simple_double( nyuko, topol, len, eff );
nlim = nyuko-1;
// fprintf( stdout, ">%s\n", name[i] );
// fprintf( stdout, "%s\n", seq[i] );
}
+ if( nunknown ) fprintf( stderr, "\nThere are %d ambiguous characters\n", nunknown );
// exit( 1 );
#if 0
pscore = 0;
for( pt=seq[i]; *pt; pt++ )
{
+// pscore += amino_dis[(int)*pt][(int)*pt];
pscore += amino_dis[(int)*pt][(int)*pt];
}
scores[i].selfscore = pscore;
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] );
+ scores[i].selfscore = localcommonsextet_p( table1, pointt[i] );
free( table1 );
}
}