source: trunk/GDE/MAFFT/mafft-7.055-with-extensions/core/contrafoldwrap.c

#include "mltaln.h"

#define DEBUG 0

static char *whereiscontrafold;

void unknown_n( char *out, char *in )
{
        while( *in )
        {
                if( *in == 'a' || *in == 'A' )
                        *out = 'A';
                else if( *in == 't' || *in == 'T' || *in == 'u' || *in == 'U' )
                        *out = 'U';
                else if( *in == 'g' || *in == 'G' )
                        *out = 'G';
                else if( *in == 'c' || *in == 'C' )
                        *out = 'C';
                else if( *in == '-' )
                        *out = '-';
                else
                        *out = 'N';

                out++;
                in++;
        }
        *out = 0;
}

void outcontrafold( FILE *fp, RNApair **pairprob, int length )
{
        int i;
        RNApair *pt;
        for( i=0; i<length; i++ ) for( pt=pairprob[i]; pt->bestpos!=-1; pt++ )
        {
                if( pt->bestpos > i ) 
                        fprintf( fp, "%d %d %f\n", i, pt->bestpos, pt->bestscore );
        }
}

#if 1
static void readcontrafold( FILE *fp, RNApair **pairprob, int length )
{
        char gett[10000];
        int *pairnum;
        char *pt;
        int i;
        int left, right;
        float prob;

        pairnum = (int *)calloc( length, sizeof( int ) );
        for( i=0; i<length; i++ ) pairnum[i] = 0;

        while( 1 )
        {
                if( feof( fp ) ) break;
                fgets( gett, 9999, fp );

//              fprintf( stderr, "gett=%s\n", gett );

                pt = gett;

                sscanf( gett, "%d ", &left );
                left--;

//              fprintf( stderr, "left=%d\n", left );
                pt = strchr( pt, ' ' ) + 1;
//              fprintf( stderr, "pt=%s\n", pt );

                while( (pt = strchr( pt, ' ' ) ) )
                {
                        pt++;
//                      fprintf( stderr, "pt=%s\n", pt );
                        sscanf( pt, "%d:%f", &right, &prob );
                        right--;

//                      fprintf( stderr, "%d-%d, %f\n", left, right, prob );

                        pairprob[left] = (RNApair *)realloc( pairprob[left], (pairnum[left]+2) * sizeof( RNApair ) );
                        pairprob[left][pairnum[left]].bestscore = prob;
                        pairprob[left][pairnum[left]].bestpos = right;
                        pairnum[left]++;
                        pairprob[left][pairnum[left]].bestscore = -1.0;
                        pairprob[left][pairnum[left]].bestpos = -1;
//                      fprintf( stderr, "%d-%d, %f\n", left, right, prob );

                        pairprob[right] = (RNApair *)realloc( pairprob[right], (pairnum[right]+2) * sizeof( RNApair ) );
                        pairprob[right][pairnum[right]].bestscore = prob;
                        pairprob[right][pairnum[right]].bestpos = left;
                        pairnum[right]++;
                        pairprob[right][pairnum[right]].bestscore = -1.0;
                        pairprob[right][pairnum[right]].bestpos = -1;
//                      fprintf( stderr, "%d-%d, %f\n", right, left, prob );
                }
        }
        free( pairnum );
}
#endif

void arguments( int argc, char *argv[] )
{
    int c;
        inputfile = NULL;
        dorp = NOTSPECIFIED;
        kimuraR = NOTSPECIFIED;
        pamN = NOTSPECIFIED;
        whereiscontrafold = NULL;

    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 'd':
                                        whereiscontrafold = *++argv;
                                        fprintf( stderr, "whereiscontrafold = %s\n", whereiscontrafold );
                                        --argc;
                                        goto nextoption;
                default:
                    fprintf( stderr, "illegal option %c\n", c );
                    argc = 0;
                    break;
            }
                }
                nextoption:
                        ;
        }
    if( argc != 0 ) 
    {
        fprintf( stderr, "options: Check source file !\n" );
        exit( 1 );
    }
}


int main( int argc, char *argv[] )
{
        static char com[10000];
        static int  *nlen;      
        int left, right;
        int res;
        static char **name, **seq, **nogap;
        static int **gapmap;
        static int *order;
        int i, j;
        FILE *infp;
        RNApair ***pairprob;
        RNApair **alnpairprob;
        RNApair *pairprobpt;
        RNApair *pt;
        int *alnpairnum;
        float prob;
        int adpos;

        arguments( argc, argv );

        if( inputfile )
        {
                infp = fopen( inputfile, "r" );
                if( !infp )
                {
                        fprintf( stderr, "Cannot open %s\n", inputfile );
                        exit( 1 );
                }
        }
        else
                infp = stdin;

        if( !whereiscontrafold )
                whereiscontrafold = "";

        getnumlen( infp );
        rewind( infp );

        if( dorp != 'd' )
        {
                fprintf( stderr, "nuc only\n" );
                exit( 1 );
        }

        seq = AllocateCharMtx( njob, nlenmax*2+1 );
        nogap = AllocateCharMtx( njob, nlenmax*2+1 );
        gapmap = AllocateIntMtx( njob, nlenmax*2+1 );
        order = AllocateIntVec( njob );
        name = AllocateCharMtx( njob, B+1 );
    nlen = AllocateIntVec( njob );
        pairprob = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
        alnpairprob = (RNApair **)calloc( nlenmax, sizeof( RNApair * ) );
        alnpairnum = AllocateIntVec( nlenmax );

        for( i=0; i<nlenmax; i++ ) alnpairnum[i] = 0;

        readData_pointer( infp, name, nlen, seq );
        fclose( infp );

        for( i=0; i<njob; i++ )
        {
                pairprob[i] = (RNApair **)calloc( nlenmax, sizeof( RNApair * ) );
                for( j=0; j<nlenmax; j++ )
                {
                        pairprob[i][j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
                        pairprob[i][j][0].bestpos = -1;
                        pairprob[i][j][0].bestscore = -1.0;
                }
                unknown_n( nogap[i], seq[i] );
                order[i] = i;
        }
        for( j=0; j<nlenmax; j++ )
        {
                alnpairprob[j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
                alnpairprob[j][0].bestpos = -1;
                alnpairprob[j][0].bestscore = -1.0;
        }


        constants( njob, seq );

        fprintf( stderr, "running contrafold\n" );
        for( i=0; i<njob; i++ )
        {
                fprintf( stderr, "%d / %d\n", i+1, njob );
                commongappick_record( 1, nogap+i, gapmap[i] );
                infp = fopen( "_contrafoldin", "w" );
                fprintf( infp, ">in\n%s\n", nogap[i] );
                fclose( infp );
#if 0 // contrafold v1
                sprintf( com, "env PATH=%s contrafold predict _contrafoldin --posteriors 0.01 > _contrafoldout", whereiscontrafold );
#else // contrafold v2
                sprintf( com, "env PATH=%s contrafold predict _contrafoldin --posteriors 0.01   _contrafoldout", whereiscontrafold );
#endif
                res = system( com );
                if( res )
                {
                        fprintf( stderr, "error in contrafold\n" );
                        fprintf( stderr, "=================================================================\n" );
                        fprintf( stderr, "=================================================================\n" );
                        fprintf( stderr, "==\n" );
                        fprintf( stderr, "== This version of MAFFT supports CONTRAfold v2.02.\n" );
                        fprintf( stderr, "== If you have a lower version of CONTRAfold installed in the\n" );
                        fprintf( stderr, "== %s directory,\n", whereiscontrafold );
                        fprintf( stderr, "== please update it!\n" );
                        fprintf( stderr, "==\n" );
                        fprintf( stderr, "=================================================================\n" );
                        fprintf( stderr, "=================================================================\n" );
                        exit( 1 );
                }


                infp = fopen( "_contrafoldout", "r" );
                readcontrafold( infp, pairprob[i], nlenmax );
                fclose( infp );
                fprintf( stdout, ">%d\n", i );
                outcontrafold( stdout, pairprob[i], nlenmax );
        }

        for( i=0; i<njob; i++ )
        {
                for( j=0; j<nlen[i]; j++ ) for( pairprobpt=pairprob[i][j]; pairprobpt->bestpos!=-1; pairprobpt++ )
                {
                        left = gapmap[i][j];
                        right = gapmap[i][pairprobpt->bestpos];
                        prob = pairprobpt->bestscore;

                        for( pt=alnpairprob[left]; pt->bestpos!=-1; pt++ )
                                if( pt->bestpos == right ) break;

                        if( pt->bestpos == -1 )
                        {
                                alnpairprob[left] = (RNApair *)realloc( alnpairprob[left], (alnpairnum[left]+2) * sizeof( RNApair ) );
                                adpos = alnpairnum[left];
                                alnpairnum[left]++;
                                alnpairprob[left][adpos].bestscore = 0.0;
                                alnpairprob[left][adpos].bestpos = right;
                                alnpairprob[left][adpos+1].bestscore = -1.0;
                                alnpairprob[left][adpos+1].bestpos = -1;
                                pt = alnpairprob[left]+adpos;
                        }
                        else
                                adpos = pt-alnpairprob[left];

                        pt->bestscore += prob;
                        if( pt->bestpos != right )
                        {
                                fprintf( stderr, "okashii!\n" );
                                exit( 1 );
                        }
//                      fprintf( stderr, "adding %d-%d, %f\n", left, right, prob );
                }
        }
        return( 0 );

#if 0
        fprintf( stdout, "result=\n" );

        for( i=0; i<nlenmax; i++ ) for( pairprobpt=alnpairprob[i]; pairprobpt->bestpos!=-1; pairprobpt++ )
        {
                pairprobpt->bestscore /= (float)njob;
                left = i;
                right = pairprobpt->bestpos;
                prob = pairprobpt->bestscore;
                fprintf( stdout, "%d-%d, %f\n", left, right, prob );
        }

        return( 0 );
#endif
}