binaries/src/mafft/core/mccaskillwrap.c

   1 #include "mltaln.h"
   2
   3 #define DEBUG 0
   4
   5 static char *whereismccaskillmea;
   6
   7 void outmccaskill( FILE *fp, RNApair **pairprob, int length )
   8 {
   9         int i;
  10         RNApair *pt;
  11         for( i=0; i<length; i++ ) for( pt=pairprob[i]; pt->bestpos!=-1; pt++ )
  12         {
  13                 if( pt->bestpos > i )
  14                         fprintf( fp, "%d %d %50.40f\n", i, pt->bestpos, pt->bestscore );
  15         }
  16 }
  17
  18 #if 1
  19 static void readrawmccaskill( FILE *fp, RNApair **pairprob, int length )
  20 {
  21         char gett[1000];
  22         int *pairnum;
  23         int i;
  24         int left, right;
  25         float prob;
  26
  27         pairnum = (int *)calloc( length, sizeof( int ) );
  28         for( i=0; i<length; i++ ) pairnum[i] = 0;
  29
  30         while( 1 )
  31         {
  32                 fgets( gett, 999, fp );
  33                 if( feof( fp ) ) break;
  34                 if( gett[0] == '>' ) continue;
  35                 sscanf( gett, "%d %d %f", &left, &right, &prob );
  36                 if( prob < 0.01 ) continue; // mxscarna to mafft ryoho ni eikyou
  37 //fprintf( stderr, "gett = %s\n", gett );
  38
  39                 if( left != right && prob > 0.0 )
  40                 {
  41                         pairprob[left] = (RNApair *)realloc( pairprob[left], (pairnum[left]+2) * sizeof( RNApair ) );
  42                         pairprob[left][pairnum[left]].bestscore = prob;
  43                         pairprob[left][pairnum[left]].bestpos = right;
  44                         pairnum[left]++;
  45                         pairprob[left][pairnum[left]].bestscore = -1.0;
  46                         pairprob[left][pairnum[left]].bestpos = -1;
  47 //                      fprintf( stderr, "%d-%d, %f\n", left, right, prob );
  48
  49                         pairprob[right] = (RNApair *)realloc( pairprob[right], (pairnum[right]+2) * sizeof( RNApair ) );
  50                         pairprob[right][pairnum[right]].bestscore = prob;
  51                         pairprob[right][pairnum[right]].bestpos = left;
  52                         pairnum[right]++;
  53                         pairprob[right][pairnum[right]].bestscore = -1.0;
  54                         pairprob[right][pairnum[right]].bestpos = -1;
  55 //                      fprintf( stderr, "%d-%d, %f\n", right, left, prob );
  56                 }
  57         }
  58         free( pairnum );
  59 }
  60 #endif
  61
  62 void arguments( int argc, char *argv[] )
  63 {
  64     int c;
  65         inputfile = NULL;
  66         dorp = NOTSPECIFIED;
  67         kimuraR = NOTSPECIFIED;
  68         pamN = NOTSPECIFIED;
  69         whereismccaskillmea = NULL;
  70
  71     while( --argc > 0 && (*++argv)[0] == '-' )
  72         {
  73         while ( (c = *++argv[0]) )
  74                 {
  75             switch( c )
  76             {
  77                                 case 'i':
  78                                         inputfile = *++argv;
  79                                         fprintf( stderr, "inputfile = %s\n", inputfile );
  80                                         --argc;
  81                                         goto nextoption;
  82                                 case 'd':
  83                                         whereismccaskillmea = *++argv;
  84                                         fprintf( stderr, "whereismccaskillmea = %s\n", whereismccaskillmea );
  85                                         --argc;
  86                                         goto nextoption;
  87                 default:
  88                     fprintf( stderr, "illegal option %c\n", c );
  89                     argc = 0;
  90                     break;
  91             }
  92                 }
  93                 nextoption:
  94                         ;
  95         }
  96     if( argc != 0 )
  97     {
  98         fprintf( stderr, "options: Check source file !\n" );
  99         exit( 1 );
 100     }
 101 }
 102
 103
 104 int main( int argc, char *argv[] )
 105 {
 106         static char com[10000];
 107         static int  *nlen;
 108         int left, right;
 109         int res;
 110         static char **name, **seq, **nogap;
 111         static int **gapmap;
 112         static int *order;
 113         int i, j;
 114         FILE *infp;
 115         RNApair ***pairprob;
 116         RNApair **alnpairprob;
 117         RNApair *pairprobpt;
 118         RNApair *pt;
 119         int *alnpairnum;
 120         float prob;
 121         int adpos;
 122
 123         arguments( argc, argv );
 124
 125         if( inputfile )
 126         {
 127                 infp = fopen( inputfile, "r" );
 128                 if( !infp )
 129                 {
 130                         fprintf( stderr, "Cannot open %s\n", inputfile );
 131                         exit( 1 );
 132                 }
 133         }
 134         else
 135                 infp = stdin;
 136
 137         if( !whereismccaskillmea )
 138                 whereismccaskillmea = "";
 139
 140         getnumlen( infp );
 141         rewind( infp );
 142
 143         if( dorp != 'd' )
 144         {
 145                 fprintf( stderr, "nuc only\n" );
 146                 exit( 1 );
 147         }
 148
 149         seq = AllocateCharMtx( njob, nlenmax*2+1 );
 150         nogap = AllocateCharMtx( njob, nlenmax*2+1 );
 151         gapmap = AllocateIntMtx( njob, nlenmax*2+1 );
 152         order = AllocateIntVec( njob );
 153         name = AllocateCharMtx( njob, B+1 );
 154     nlen = AllocateIntVec( njob );
 155         pairprob = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
 156         alnpairprob = (RNApair **)calloc( nlenmax, sizeof( RNApair * ) );
 157         alnpairnum = AllocateIntVec( nlenmax );
 158
 159         for( i=0; i<nlenmax; i++ ) alnpairnum[i] = 0;
 160
 161         readData_pointer( infp, name, nlen, seq );
 162         fclose( infp );
 163
 164         for( i=0; i<njob; i++ )
 165         {
 166                 pairprob[i] = (RNApair **)calloc( nlenmax, sizeof( RNApair * ) );
 167                 for( j=0; j<nlenmax; j++ )
 168                 {
 169                         pairprob[i][j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
 170                         pairprob[i][j][0].bestpos = -1;
 171                         pairprob[i][j][0].bestscore = -1.0;
 172                 }
 173                 strcpy( nogap[i], seq[i] );
 174                 order[i] = i;
 175         }
 176         for( j=0; j<nlenmax; j++ )
 177         {
 178                 alnpairprob[j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
 179                 alnpairprob[j][0].bestpos = -1;
 180                 alnpairprob[j][0].bestscore = -1.0;
 181         }
 182
 183
 184         constants( njob, seq );
 185
 186         fprintf( stderr, "Running mxscarna with the mccaskill_mea mode.\n" );
 187         for( i=0; i<njob; i++ )
 188         {
 189                 fprintf( stderr, "%d / %d\n", i+1, njob );
 190                 commongappick_record( 1, nogap+i, gapmap[i] );
 191                 infp = fopen( "_mccaskillinorg", "w" );
 192 //              fprintf( infp, ">in\n%s\n", nogap[i] );
 193                 fprintf( infp, ">in\n" );
 194                 write1seq( infp, nogap[i] );
 195                 fclose( infp );
 196
 197                 system( "tr -d '\\r' < _mccaskillinorg > _mccaskillin" ); // for cygwin, wakaran
 198                 sprintf( com, "env PATH=%s mxscarnamod -m -writebpp  _mccaskillin > _mccaskillout 2>_dum", whereismccaskillmea );
 199                 res = system( com );
 200
 201                 if( res )
 202                 {
 203                         fprintf( stderr, "ERROR IN mccaskill_mea\n" );
 204                         exit( 1 );
 205                 }
 206
 207                 infp = fopen( "_mccaskillout", "r" );
 208                 readrawmccaskill( infp, pairprob[i], nlenmax );
 209                 fclose( infp );
 210                 fprintf( stdout, ">%d\n", i );
 211                 outmccaskill( stdout, pairprob[i], nlenmax );
 212         }
 213
 214         for( i=0; i<njob; i++ )
 215         {
 216                 for( j=0; j<nlen[i]; j++ ) for( pairprobpt=pairprob[i][j]; pairprobpt->bestpos!=-1; pairprobpt++ )
 217                 {
 218                         left = gapmap[i][j];
 219                         right = gapmap[i][pairprobpt->bestpos];
 220                         prob = pairprobpt->bestscore;
 221
 222                         for( pt=alnpairprob[left]; pt->bestpos!=-1; pt++ )
 223                                 if( pt->bestpos == right ) break;
 224
 225                         if( pt->bestpos == -1 )
 226                         {
 227                                 alnpairprob[left] = (RNApair *)realloc( alnpairprob[left], (alnpairnum[left]+2) * sizeof( RNApair ) );
 228                                 adpos = alnpairnum[left];
 229                                 alnpairnum[left]++;
 230                                 alnpairprob[left][adpos].bestscore = 0.0;
 231                                 alnpairprob[left][adpos].bestpos = right;
 232                                 alnpairprob[left][adpos+1].bestscore = -1.0;
 233                                 alnpairprob[left][adpos+1].bestpos = -1;
 234                                 pt = alnpairprob[left]+adpos;
 235                         }
 236                         else
 237                                 adpos = pt-alnpairprob[left];
 238
 239                         pt->bestscore += prob;
 240                         if( pt->bestpos != right )
 241                         {
 242                                 fprintf( stderr, "okashii!\n" );
 243                                 exit( 1 );
 244                         }
 245 //                      fprintf( stderr, "adding %d-%d, %f\n", left, right, prob );
 246                 }
 247         }
 248
 249         for( i=0; i<njob; i++ )
 250         {
 251                 for( j=0; j<nlenmax; j++ ) free( pairprob[i][j] );
 252                 free( pairprob[i] );
 253         }
 254         free( pairprob );
 255         for( j=0; j<nlenmax; j++ ) free( alnpairprob[j] );
 256         free( alnpairprob );
 257         free( alnpairnum );
 258         return( 0 );
 259
 260 #if 0
 261         fprintf( stdout, "result=\n" );
 262
 263         for( i=0; i<nlenmax; i++ ) for( pairprobpt=alnpairprob[i]; pairprobpt->bestpos!=-1; pairprobpt++ )
 264         {
 265                 pairprobpt->bestscore /= (float)njob;
 266                 left = i;
 267                 right = pairprobpt->bestpos;
 268                 prob = pairprobpt->bestscore;
 269                 fprintf( stdout, "%d-%d, %f\n", left, right, prob );
 270         }
 271
 272         return( 0 );
 273 #endif
 274 }