JWS-112 Bumping version of Mafft to version 7.310.
[jabaws.git] / binaries / src / mafft / core / pairash.c
index 1e83e72..a18267c 100644 (file)
@@ -6,6 +6,8 @@
 
 static int usecache;
 static char *whereispairalign;
+static char *odir;
+static char *pdir;
 static double scale;
 static int *alreadyoutput;
 static int equivthreshold;
@@ -327,12 +329,12 @@ static int checkcbeta( FILE *fp )
 }
 
 
-static float calltmalign( char **mseq1, char **mseq2, double *equiv, char *fname1, char *chain1, char *fname2, char *chain2, int alloclen )
+static double calltmalign( char **mseq1, char **mseq2, double *equiv, char *fname1, char *chain1, char *fname2, char *chain2, int alloclen )
 {
        FILE *fp;
        int res;
        static char com[10000];
-       float value;
+       double value;
        char cachedir[10000];
        char cachefile[10000];
        int runnow;
@@ -399,17 +401,17 @@ static float calltmalign( char **mseq1, char **mseq2, double *equiv, char *fname
 //     fprintf( stderr, "*mseq1 = %s\n", *mseq1 );
 //     fprintf( stderr, "*mseq2 = %s\n", *mseq2 );
 
-       value = (float)naivepairscore11( *mseq1, *mseq2, penalty );
+       value = (double)naivepairscore11( *mseq1, *mseq2, penalty );
 
        return( value );
 }
 
-static float callrash( int mem1, int mem2, char **mseq1, char **mseq2, double *equiv, char *fname1, char *chain1, char *fname2, char *chain2, int alloclen )
+static double callrash( int mem1, int mem2, char **mseq1, char **mseq2, double *equiv, char *fname1, char *fname2, int alloclen )
 {
        FILE *fp;
-       int res;
+//     int res;
        static char com[10000];
-       float value;
+       double value;
        char cachedir[10000];
        char cachefile[10000];
        int runnow;
@@ -420,28 +422,31 @@ static float callrash( int mem1, int mem2, char **mseq1, char **mseq2, double *e
 
        if( usecache )
        {
-               sprintf( cachedir, "%s/.rashoutcache", getenv( "HOME" ) );
-               sprintf( com, "mkdir -p %s", cachedir );
-               system( com );
+//             sprintf( cachedir, "tmp" );
+               sprintf( cachedir, "%s", pdir );
 
-               sprintf( cachefile, "%s/%s%s-%s%s", cachedir, fname1, chain1, fname2, chain2 );
+               sprintf( cachefile, "%s/%s.%s.rash", cachedir, fname1, fname2 );
 
                runnow = 0;
                fp = fopen( cachefile, "r" );
-               if( fp == NULL ) runnow = 1;
+               if( fp == NULL )
+               {
+                       fprintf( stderr, "Cannot open %s\n", cachefile );
+                       exit( 1 );
+               }
                else
                {
-                       fgets( com, 100, fp );
-                       if( strncmp( com, "successful", 10 ) ) runnow = 1;
                        fclose( fp );
                }
        }
        else
        {
-               runnow = 1;
+               fprintf( stderr, "Not supported!\n" );
+               exit( 1 );
        }
 
-       if( runnow )
+#if 0
+       if( 0 )
        {
 #if 0
                sprintf( com, "ln -s %s %s.pdb 2>_dum", fname1, fname1 );
@@ -473,18 +478,17 @@ static float callrash( int mem1, int mem2, char **mseq1, char **mseq2, double *e
        
        }
        else
+#endif
        {
-               fprintf( stderr, "Use cache!\n" );
-               sprintf( com, "grep -v successful %s > %s.rashout", cachefile, pairid );
+               fprintf( stderr, "Use cache! cachefile = %s\n", cachefile );
+               sprintf( com, "cat %s > %s.rashout", cachefile, pairid );
                system( com );
        }
 
        if( usecache && runnow )
        {
-               sprintf( com, "echo successful >  %s", cachefile );
-               system( com );
-               sprintf( com, "cat %s.rashout >>  %s", pairid, cachefile );
-               system( com );
+               fprintf( stderr, "Okashii! usechache=%d, runnow=%d\n", usecache, runnow );
+               exit( 1 );
        }
 
        sprintf( com, "%s.rashout", pairid );
@@ -503,7 +507,7 @@ static float callrash( int mem1, int mem2, char **mseq1, char **mseq2, double *e
 //     fprintf( stderr, "*mseq2 = %s\n", *mseq2 );
 
 
-       value = (float)naivepairscore11( *mseq1, *mseq2, penalty );
+       value = (double)naivepairscore11( *mseq1, *mseq2, penalty );
 
        return( value );
 }
@@ -518,7 +522,7 @@ static void preparetmalign( FILE *fp, char ***strfiles, char ***chainids, char *
        int linenum, istr, nstr;
        FILE *checkfp;
        char *sline; 
-       int use[1000];
+       int use[10000];
        linenum = 0;
        nstr = 0;
        while( 1 )
@@ -615,7 +619,7 @@ static void preparetmalign( FILE *fp, char ***strfiles, char ***chainids, char *
        }
 }
 
-static void prepareash( FILE *fp, char ***strfiles, char ***chainids, char ***seqpt, char ***mseq1pt, char ***mseq2pt, double **equivpt, int *alloclenpt )
+static void prepareash( FILE *fp, char *inputfile, char ***strfiles, char ***chainids, char ***seqpt, char ***mseq1pt, char ***mseq2pt, double **equivpt, int *alloclenpt )
 {
        int i, res;
        char *dumseq;
@@ -623,11 +627,13 @@ static void prepareash( FILE *fp, char ***strfiles, char ***chainids, char ***se
        char fname[1000];
        char command[1000];
        int linenum, istr, nstr;
-       FILE *checkfp;
+//     FILE *checkfp;
        char *sline; 
-       int use[1000];
+       int use[10000];
        linenum = 0;
        nstr = 0;
+
+       fprintf( stderr, "inputfile = %s\n", inputfile );
        while( 1 )
        {
                fgets( line, 999, fp );
@@ -641,20 +647,13 @@ static void prepareash( FILE *fp, char ***strfiles, char ***chainids, char ***se
                        continue;
                }
                extractfirstword( sline );
+#if 0
                checkfp = fopen( sline, "r" );
                if( checkfp == NULL )
                {
                        fprintf( stderr, "Cannot open %s.\n", sline );
                        exit( 1 );
                }
-#if 0
-               fgets( linec, 999, checkfp );
-               if( strncmp( "HEADER ", linec, 7 ) )
-               {
-                       fprintf( stderr, "Check the format of %s.\n", sline );
-                       exit( 1 );
-               }
-#endif
                if( checkcbeta( checkfp ) ) 
                {
                        fprintf( stderr, "%s has no C-beta atoms.\n", sline );
@@ -663,6 +662,9 @@ static void prepareash( FILE *fp, char ***strfiles, char ***chainids, char ***se
                else
                        nstr++;
                fclose( checkfp );
+#else
+               nstr++;
+#endif
                linenum++;
        }
        njob = nstr;
@@ -679,6 +681,7 @@ static void prepareash( FILE *fp, char ***strfiles, char ***chainids, char ***se
                fgets( line, 999, fp );
                if( feof( fp ) ) break;
                sline = strip( line );
+               fprintf( stderr, "sline = %s\n", sline );
                if( use[linenum++] ) 
                {
                        (*chainids)[istr][0] = getchainid( sline );
@@ -686,6 +689,7 @@ static void prepareash( FILE *fp, char ***strfiles, char ***chainids, char ***se
                        extractfirstword( sline );
                        sprintf( fname, "%s", sline );
                        cutpath( fname );
+#if 0
                        sprintf( command, "cp %s %s.pdb", sline, fname );
                        system( command );
                        sprintf( command, "perl \"%s/clean.pl\" %s.pdb", whereispairalign, fname );
@@ -695,6 +699,7 @@ static void prepareash( FILE *fp, char ***strfiles, char ***chainids, char ***se
                                fprintf( stderr, "error: Install clean.pl\n" );
                                exit( 1 );
                        }
+#endif
                        strcpy( (*strfiles)[istr++], fname );
                }
        }
@@ -708,6 +713,18 @@ static void prepareash( FILE *fp, char ***strfiles, char ***chainids, char ***se
        alreadyoutput = AllocateIntVec( njob );
        for( i=0; i<njob; i++ ) alreadyoutput[i] = 0;
 
+       fprintf( stderr, "Running pdp_ash_batch.pl..\n" );
+//     sprintf( command, "/opt/protein/share/domains/code/pdp_ash/pdp_ash_batch.pl -f %s -d tmp -i %d", inputfile, wheretooutput  );
+       sprintf( command, "/opt/protein/share/mafftash/pdp_ash/pdp_ash_batch.pl -f %s -d %s -i %s", inputfile, pdir, odir );
+       res = system( command );
+       if( res )
+       {
+               fprintf( stderr, "Ask KM!\n" );
+               exit( 1 );
+       }
+       fprintf( stderr, "done\n" );
+
+
        for( i=0; i<istr; i++ )
        {
                fprintf( stderr, "i=%d\n", i );
@@ -716,8 +733,8 @@ static void prepareash( FILE *fp, char ***strfiles, char ***chainids, char ***se
                (*mseq1pt)[0] = (*seqpt)[i];
                (*mseq2pt)[0] = dumseq;
 
-               callrash( i, i, *mseq1pt, *mseq2pt, *equivpt, (*strfiles)[i], (*chainids)[i], (*strfiles)[i], (*chainids)[i], *alloclenpt );
-               fprintf( stdout, ">%d_%s-%s\n%s\n", i+1, (*strfiles)[i], (*chainids)[i], (*seqpt)[i] );
+               callrash( i, i, *mseq1pt, *mseq2pt, *equivpt, (*strfiles)[i], (*strfiles)[i], *alloclenpt );
+               fprintf( stdout, ">%d_%s\n%s\n", i+1, (*strfiles)[i], (*seqpt)[i] );
                alreadyoutput[i] = 1;
        }
 }
@@ -776,6 +793,8 @@ void arguments( int argc, char *argv[] )
        fftThreshold = NOTSPECIFIED;
        RNAppenalty = NOTSPECIFIED;
        RNApthr = NOTSPECIFIED;
+       odir = "";
+       pdir = "";
 
     while( --argc > 0 && (*++argv)[0] == '-' )
        {
@@ -809,25 +828,25 @@ void arguments( int argc, char *argv[] )
                                        --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 );
@@ -838,16 +857,26 @@ void arguments( int argc, char *argv[] )
                                        fprintf( stderr, "whereispairalign = %s\n", whereispairalign );
                                        --argc; 
                                        goto nextoption;
+                               case 'o':
+                                       odir = *++argv;
+                                       fprintf( stderr, "odir = %s\n", odir );
+                                       --argc; 
+                                       goto nextoption;
+                               case 'p':
+                                       pdir = *++argv;
+                                       fprintf( stderr, "pdir = %s\n", pdir );
+                                       --argc; 
+                                       goto nextoption;
                                case 't':
-                                       equivthreshold = atoi( *++argv );
+                                       equivthreshold = myatoi( *++argv );
                                        --argc;
                                        goto nextoption;
                                case 'w':
-                                       equivwinsize = atoi( *++argv );
+                                       equivwinsize = myatoi( *++argv );
                                        --argc;
                                        goto nextoption;
                                case 'l':
-                                       equivshortestlen = atoi( *++argv );
+                                       equivshortestlen = myatoi( *++argv );
                                        --argc;
                                        goto nextoption;
                                case 's':
@@ -942,11 +971,11 @@ void arguments( int argc, char *argv[] )
 /* modification end. */
 #if 0
                                case 'z':
-                                       fftThreshold = atoi( *++argv );
+                                       fftThreshold = myatoi( *++argv );
                                        --argc; 
                                        goto nextoption;
                                case 'w':
-                                       fftWinSize = atoi( *++argv );
+                                       fftWinSize = myatoi( *++argv );
                                        --argc;
                                        goto nextoption;
                                case 'Z':
@@ -992,12 +1021,12 @@ int countamino( char *s, int end )
        return( val );
 }
 
-static void pairalign( char name[M][B], int nlen[M], char **seq, char **aseq, char **mseq1, char **mseq2, double *equiv, double *effarr, char **strfiles, char **chainids, int alloclen )
+static void pairalign( char **name, int nlen[M], char **seq, char **aseq, char **mseq1, char **mseq2, double *equiv, double *effarr, char **strfiles, char **chainids, int alloclen )
 {
        int i, j, ilim;
        int clus1, clus2;
        int off1, off2;
-       float pscore = 0.0; // by D.Mathog
+       double pscore = 0.0; // by D.Mathog
        static char *indication1, *indication2;
        FILE *hat2p, *hat3p;
        static double **distancemtx;
@@ -1085,8 +1114,8 @@ static void pairalign( char name[M][B], int nlen[M], char **seq, char **aseq, ch
 
                        strcpy( aseq[i], seq[i] );
                        strcpy( aseq[j], seq[j] );
-                       clus1 = conjuctionfortbfast( pair, i, aseq, mseq1, effarr1, effarr, indication1 );
-                       clus2 = conjuctionfortbfast( pair, j, aseq, mseq2, effarr2, effarr, indication2 );
+                       clus1 = conjuctionfortbfast_old( pair, i, aseq, mseq1, effarr1, effarr, indication1 );
+                       clus2 = conjuctionfortbfast_old( pair, j, aseq, mseq2, effarr2, effarr, indication2 );
        //              fprintf( stderr, "mseq1 = %s\n", mseq1[0] );
        //              fprintf( stderr, "mseq2 = %s\n", mseq2[0] );
        
@@ -1109,7 +1138,7 @@ static void pairalign( char name[M][B], int nlen[M], char **seq, char **aseq, ch
                                                break;
                                        case( 'R' ):
                                                fprintf( stderr, "  Calling PDP_ASH.pl %d-%d/%d    \r", i+1, j+1, njob );
-                                               pscore = callrash( i, j, mseq1, mseq2, equiv, strfiles[i], chainids[i], strfiles[j], chainids[j], alloclen );
+                                               pscore = callrash( i, j, mseq1, mseq2, equiv, strfiles[i], strfiles[j], alloclen );
                                                off1 = off2 = 0;
                                                break;
                                        ErrorExit( "ERROR IN SOURCE FILE" );
@@ -1121,14 +1150,14 @@ static void pairalign( char name[M][B], int nlen[M], char **seq, char **aseq, ch
                        fprintf( stderr, "score = %10.2f (%d,%d)\n", pscore, i, j );
 #endif
 
-                       putlocalhom_str( mseq1[0], mseq2[0], equiv, scale, localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
+                       putlocalhom_str( mseq1[0], mseq2[0], equiv, scale, localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ), 'k' );
 #if 1
 
                        if( alreadyoutput[i] == 0 )
                        {
                                alreadyoutput[i] = 1;
                                gappick0( seq[i], mseq1[0] );
-                               fprintf( stdout, ">%d_%s-%s\n%s\n", i+1, strfiles[i], chainids[i], seq[i] );
+                               fprintf( stdout, ">%d_%s\n%s\n", i+1, strfiles[i], seq[i] );
                                strcpy( checkseq[i], seq[i] );
                        }
                        else
@@ -1195,7 +1224,7 @@ static void pairalign( char name[M][B], int nlen[M], char **seq, char **aseq, ch
 
        hat2p = fopen( hat2file, "w" );
        if( !hat2p ) ErrorExit( "Cannot open hat2." );
-       WriteHat2( hat2p, njob, name, distancemtx );
+       WriteHat2_pointer( hat2p, njob, name, distancemtx );
        fclose( hat2p );
 
        fprintf( stderr, "##### writing hat3\n" );
@@ -1285,7 +1314,7 @@ static void WriteOptions( FILE *fp )
 int main( int argc, char *argv[] )
 {
        static int  nlen[M];    
-       static char name[M][B], **seq;
+       static char **name, **seq;
        static char **mseq1, **mseq2;
        static char **aseq;
        static char **bseq;
@@ -1328,12 +1357,13 @@ int main( int argc, char *argv[] )
        nlenmax = 10000; // tekitou
 
        if( alg == 'R' )
-               prepareash( infp, &strfiles, &chainids, &seq, &mseq1, &mseq2, &equiv, &alloclen );
+               prepareash( infp, inputfile, &strfiles, &chainids, &seq, &mseq1, &mseq2, &equiv, &alloclen );
        else if( alg == 'T' )
                preparetmalign( infp, &strfiles, &chainids, &seq, &mseq1, &mseq2, &equiv, &alloclen );
 
        fclose( infp );
 
+       name = AllocateCharMtx( njob, B+1 );
        aseq = AllocateCharMtx( njob, nlenmax*2+1 );
        bseq = AllocateCharMtx( njob, nlenmax*2+1 );
        eff = AllocateDoubleVec( njob );