website/archive/binaries/mac/src/mafft/core/addfunctions.c

   1 #include "mltaln.h"
   2
   3 void profilealignment( int n0, int n1, int n2, char **aln0, char **aln1, char **aln2, int alloclen, char alg ) // n1 ha allgap
   4 {
   5         int i, newlen;
   6         double *effarr0, *effarr2;
   7         float dumfl;
   8         double eff;
   9         effarr0 = AllocateDoubleVec( n0 );
  10         effarr2 = AllocateDoubleVec( n2 );
  11
  12         commongappick( n0, aln0 );
  13         commongappick( n2, aln2 );
  14
  15
  16         eff = 1.0 / (double)n0; for( i=0; i<n0; i++ ) effarr0[i] = eff;
  17         eff = 1.0 / (double)n2; for( i=0; i<n2; i++ ) effarr2[i] = eff;
  18
  19         if( alg == 'M' )
  20                 MSalignmm( aln0, aln2, effarr0, effarr2, n0, n2, alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 ); //outgap=1??
  21         else
  22                 A__align( aln0, aln2, effarr0, effarr2, n0, n2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 ); //outgap=1??
  23
  24         newlen = strlen( aln0[0] );
  25
  26         for( i=0; i<newlen; i++ ) aln1[0][i] = '-';
  27         aln1[0][i] = 0;
  28         for( i=1; i<n1; i++ ) strcpy( aln1[i], aln1[0] );
  29
  30
  31         free( effarr0 );
  32         free( effarr2 );
  33 }
  34
  35 void eq2dash( char *s )
  36 {
  37         while( *s )
  38         {
  39                 if( *s == '=' ) *s = '-';
  40                 s++;
  41         }
  42 }
  43
  44 void findnewgaps( int n, char **seq, int *gaplen )
  45 {
  46         int i, pos, len;
  47         len = strlen( seq[0] );
  48
  49 //      fprintf( stderr, "seq[0] = %s\n", seq[0] );
  50         for( i=0; i<len; i++ ) gaplen[i] = 0;
  51
  52         pos = 0;
  53         for( i=0; i<len; i++ )
  54         {
  55                 if( seq[0][i] == '=' )
  56                 {
  57 //                      fprintf( stderr, "Newgap! pos = %d\n", pos );
  58                         gaplen[pos]++;
  59                 }
  60                 else
  61                         pos++;
  62         }
  63 }
  64
  65 void findcommongaps( int n, char **seq, int *gaplen )
  66 {
  67         int i, j, pos, len, len1;
  68         len = strlen( seq[0] );
  69         len1 = len+1;
  70
  71 //      fprintf( stderr, "seq[0] = %s\n", seq[0] );
  72         for( i=0; i<len1; i++ ) gaplen[i] = 0;
  73
  74         pos = 0;
  75         for( i=0; i<len; i++ )
  76         {
  77                 for( j=0; j<n; j++ )
  78                         if( seq[j][i] != '-' ) break;
  79
  80                 if( j == n ) gaplen[pos]++;
  81                 else
  82                         pos++;
  83         }
  84 #if 0
  85         for( i=0; i<pos; i++ )
  86         {
  87                 fprintf( stderr, "vec[%d] = %d\n", i, gaplen[i] );
  88         }
  89 #endif
  90 }
  91
  92 void adjustgapmap( int newlen, int *gapmap, char *seq )
  93 {
  94         int j;
  95         int pos;
  96         int newlen1 = newlen+1;
  97         int *tmpmap;
  98
  99         tmpmap = AllocateIntVec( newlen+2 );
 100         j = 0;
 101         pos = 0;
 102         while( *seq )
 103         {
 104 //              fprintf( stderr, "j=%d *seq = %c\n", j, *seq );
 105                 if( *seq++ == '=' )
 106                         tmpmap[j++] = 0;
 107                 else
 108                 {
 109                         tmpmap[j++] = gapmap[pos++];
 110                 }
 111         }
 112         tmpmap[j++] = gapmap[pos];
 113
 114         for(j=0; j<newlen1; j++)
 115                 gapmap[j] = tmpmap[j];
 116
 117         free( tmpmap );
 118 }
 119
 120 void insertnewgaps( int njob, int *alreadyaligned, char **seq, int *ex1, int *ex2, int *gaplen, int *gapmap, int alloclen, char alg )
 121 {
 122         int *mar;
 123         char *gaps;
 124         char *cptr;
 125         int i, j, k, len, rep, len0;
 126         char **mseq2, **mseq0, **mseq1;
 127         char **aseq;
 128         int ngroup2, ngroup0, ngroup1;
 129         int *list0, *list1, *list2;
 130         int posin12, gapshift, newpos;
 131
 132         mar = calloc( njob, sizeof( int ) );
 133         list0 = calloc( njob, sizeof( int ) );
 134         list1 = calloc( njob, sizeof( int ) );
 135         list2 = calloc( njob, sizeof( int ) );
 136
 137         for( i=0; i<njob; i++ ) mar[i] = 0;
 138         for( i=0; i<njob; i++ )
 139         {
 140                 if( alreadyaligned[i]==0 ) mar[i] = 3;
 141         }
 142         for( i=0; (k=ex1[i])>-1; i++ )
 143         {
 144                 mar[k] = 1;
 145 //              fprintf( stderr, "excluding %d\n", ex1[i] );
 146         }
 147         for( i=0; (k=ex2[i])>-1; i++ )
 148         {
 149                 mar[k] = 2;
 150 //              fprintf( stderr, "excluding %d\n", ex2[i] );
 151         }
 152
 153         ngroup2 = ngroup1 = ngroup0 = 0;
 154         for( i=0; i<njob; i++ )
 155         {
 156                 if( mar[i] == 2 )
 157                 {
 158                         list2[ngroup2] = i;
 159                         ngroup2++;
 160                 }
 161                 if( mar[i] == 1 )
 162                 {
 163                         list1[ngroup1] = i;
 164                         ngroup1++;
 165                 }
 166                 if( mar[i] == 0 )
 167                 {
 168                         list0[ngroup0] = i;
 169                         ngroup0++;
 170                 }
 171         }
 172         list0[ngroup0] = list1[ngroup1] = list2[ngroup2] = -1;
 173         if( ngroup0 == 0 )
 174         {
 175                 fprintf( stderr, "Nothing to do\n" );
 176                 free( mar );
 177                 free( list0 );
 178                 free( list1 );
 179                 free( list2 );
 180                 return;
 181         }
 182
 183         for( i=0; i<njob; i++ ) if( mar[i] == 0 ) break;
 184         rep = i;
 185         len = strlen( seq[rep] );
 186         len0 = len+1;
 187
 188 //
 189 //      if( i == njob )
 190 //      {
 191 ////            fprintf( stderr, "Nothing to do\n" );
 192 //              free( mar );
 193 //              return;
 194 //      }
 195
 196         mseq2 = AllocateCharMtx( ngroup2, alloclen );
 197         mseq1 = AllocateCharMtx( ngroup1, alloclen );
 198         mseq0 = AllocateCharMtx( ngroup0, alloclen );
 199         aseq = AllocateCharMtx( njob, alloclen );
 200         gaps = calloc( alloclen, sizeof( char ) );
 201
 202         for( i=0; i<njob; i++ ) aseq[i][0] = 0;
 203         posin12 = 0;
 204         for( j=0; j<len0; j++ )
 205         {
 206                 if( gaplen[j] )
 207                 {
 208                         for( i=0; i<ngroup0; i++ ) mseq0[i][0] = 0;
 209                         for( i=0; i<ngroup1; i++ ) mseq1[i][0] = 0;
 210                         for( i=0; i<ngroup2; i++ ) mseq2[i][0] = 0;
 211
 212                         gapshift = gaplen[j];
 213                         cptr = gaps;
 214                         while( gapshift-- ) *cptr++ = '-';
 215                         *cptr = 0;
 216                         gapshift = gaplen[j];
 217
 218                         for( i=0; i<ngroup0; i++ ) strncat( mseq0[i], gaps, gapshift );
 219                         for( i=0; i<ngroup1; i++ ) strncat( mseq1[i], seq[list1[i]]+posin12, gapshift );
 220                         for( i=0; i<ngroup2; i++ ) strncat( mseq2[i], seq[list2[i]]+posin12, gapshift );
 221                         posin12 += gapshift;
 222
 223                         gapshift = gapmap[posin12];
 224 //                      fprintf( stderr, "gapmap[%d] kouho = %d\n", posin12, gapmap[posin12] );
 225
 226
 227                         for( i=0; i<ngroup0; i++ ) strncat( mseq0[i], seq[list0[i]]+j, gapshift );
 228                         for( i=0; i<ngroup1; i++ ) strncat( mseq1[i], seq[list1[i]]+posin12, gapshift );
 229                         for( i=0; i<ngroup2; i++ ) strncat( mseq2[i], seq[list2[i]]+posin12, gapshift );
 230 #if 0
 231                         for( i=0; i<ngroup0; i++ ) fprintf( stderr, "### mseq0[%d] = %s\n", i, mseq0[i] );
 232                         for( i=0; i<ngroup1; i++ ) fprintf( stderr, "### mseq1[%d] = %s\n", i, mseq1[i] );
 233                         for( i=0; i<ngroup2; i++ ) fprintf( stderr, "### mseq2[%d] = %s\n", i, mseq2[i] );
 234 #endif
 235
 236                         if( gapshift ) profilealignment( ngroup0, ngroup1, ngroup2, mseq0, mseq1, mseq2, alloclen, alg );
 237
 238                         j += gapshift;
 239                         posin12 += gapshift;
 240
 241                         for( i=0; i<ngroup0; i++ ) strcat( aseq[list0[i]], mseq0[i] );
 242                         for( i=0; i<ngroup1; i++ ) strcat( aseq[list1[i]], mseq1[i] );
 243                         for( i=0; i<ngroup2; i++ ) strcat( aseq[list2[i]], mseq2[i] );
 244                 }
 245
 246                 newpos = strlen( aseq[rep] );
 247                 for( i=0; i<ngroup0; i++ ) aseq[list0[i]][newpos] = seq[list0[i]][j];
 248                 for( i=0; i<ngroup1; i++ ) aseq[list1[i]][newpos] = seq[list1[i]][posin12];
 249                 for( i=0; i<ngroup2; i++ ) aseq[list2[i]][newpos] = seq[list2[i]][posin12];
 250                 for( i=0; i<ngroup0; i++ ) aseq[list0[i]][newpos+1] = 0;
 251                 for( i=0; i<ngroup1; i++ ) aseq[list1[i]][newpos+1] = 0;
 252                 for( i=0; i<ngroup2; i++ ) aseq[list2[i]][newpos+1] = 0;
 253
 254                 posin12++;
 255         }
 256
 257 //      for( i=0; i<njob; i++ ) if( mar[i] != 3 ) strcpy( seq[i], aseq[i] );
 258         for( i=0; i<ngroup0; i++ ) strcpy( seq[list0[i]], aseq[list0[i]] );
 259         for( i=0; i<ngroup1; i++ ) strcpy( seq[list1[i]], aseq[list1[i]] );
 260         for( i=0; i<ngroup2; i++ ) strcpy( seq[list2[i]], aseq[list2[i]] );
 261
 262         free( mar );
 263         free( gaps );
 264         free( list0 );
 265         free( list1 );
 266         free( list2 );
 267         FreeCharMtx( mseq2 );
 268         FreeCharMtx( mseq0 );
 269 }
 270
 271
 272 void restorecommongaps( int njob, char **seq, int *ex1, int *ex2, int *gaplen, int alloclen )
 273 {
 274         int *mar;
 275         char *tmpseq;
 276         char *cptr;
 277         int *iptr;
 278         int *tmpgaplen;
 279         int i, j, k, len, rep, len1;
 280
 281         mar = calloc( njob, sizeof( int ) );
 282         tmpseq = calloc( alloclen, sizeof( char ) );
 283         tmpgaplen = calloc( alloclen, sizeof( int ) );
 284 //      tmpseq = calloc( alloclen+2, sizeof( char ) );
 285 //      tmpgaplen = calloc( alloclen+2, sizeof( int ) );
 286
 287
 288         for( i=0; i<njob; i++ ) mar[i] = 0;
 289         for( i=0; (k=ex1[i])>-1; i++ )
 290         {
 291                 mar[k] = 1;
 292 //              fprintf( stderr, "excluding %d\n", ex1[i] );
 293         }
 294         for( i=0; (k=ex2[i])>-1; i++ )
 295         {
 296                 mar[k] = 1;
 297 //              fprintf( stderr, "excluding %d\n", ex2[i] );
 298         }
 299
 300         for( i=0; i<njob; i++ )
 301                 if( mar[i] ) break;
 302
 303         if( i == njob )
 304         {
 305 //              fprintf( stderr, "Nothing to do\n" );
 306                 free( mar );
 307                 free( tmpseq );
 308                 free( tmpgaplen );
 309                 return;
 310         }
 311         rep = i;
 312         len = strlen( seq[rep] );
 313         len1 = len+1;
 314
 315
 316         for( i=0; i<njob; i++ )
 317         {
 318                 if( mar[i] == 0 ) continue;
 319                 cptr = tmpseq;
 320                 for( j=0; j<len1; j++ )
 321                 {
 322                         for( k=0; k<gaplen[j]; k++ )
 323                                 *(cptr++) = '-';
 324                         *(cptr++) = seq[i][j];
 325                 }
 326                 *cptr = 0;
 327                 strcpy( seq[i], tmpseq );
 328         }
 329
 330         iptr = tmpgaplen;
 331         for( j=0; j<len1; j++ )
 332         {
 333                 *(iptr++) = gaplen[j];
 334                 for( k=0; k<gaplen[j]; k++ )
 335                         *(iptr++) = 0;
 336         }
 337         *iptr = -1;
 338
 339         iptr = tmpgaplen;
 340         while( *iptr != -1 ) *gaplen++ = *iptr++;
 341
 342         free( mar );
 343         free( tmpseq );
 344         free( tmpgaplen );
 345 }
 346
 347 #if 0
 348 int samemember( int *mem, int *cand )
 349 {
 350         int i, j;
 351
 352 #if 0
 353         fprintf( stderr, "mem = " );
 354         for( i=0; mem[i]>-1; i++ )      fprintf( stderr, "%d ", mem[i] );
 355         fprintf( stderr, "\n" );
 356
 357         fprintf( stderr, "cand = " );
 358         for( i=0; cand[i]>-1; i++ )     fprintf( stderr, "%d ", cand[i] );
 359         fprintf( stderr, "\n" );
 360 #endif
 361
 362         for( i=0, j=0; mem[i]>-1; )
 363         {
 364                 if( mem[i++] != cand[j++] ) return( 0 );
 365         }
 366
 367         if( cand[j] == -1 )
 368         {
 369                 return( 1 );
 370         }
 371         else
 372         {
 373                 return( 0 );
 374         }
 375 }
 376 #else
 377 int samemember( int *mem, int *cand )
 378 {
 379         int i, j;
 380         int nm, nc;
 381 #if 0
 382         fprintf( stderr, "mem = " );
 383         for( i=0; mem[i]>-1; i++ )      fprintf( stderr, "%d ", mem[i] );
 384         fprintf( stderr, "\n" );
 385
 386         fprintf( stderr, "cand = " );
 387         for( i=0; cand[i]>-1; i++ )     fprintf( stderr, "%d ", cand[i] );
 388         fprintf( stderr, "\n" );
 389 #endif
 390
 391         nm = 0; for( i=0; mem[i]>-1; i++ ) nm++;
 392         nc = 0; for( i=0; cand[i]>-1; i++ ) nc++;
 393
 394         if( nm != nc ) return( 0 );
 395
 396         for( i=0; mem[i]>-1; i++ )
 397         {
 398                 for( j=0; cand[j]>-1; j++ )
 399                         if( mem[i] == cand[j] ) break;
 400                 if( cand[j] == -1 ) return( 0 );
 401         }
 402
 403         if( mem[i] == -1 )
 404         {
 405                 return( 1 );
 406         }
 407         else
 408         {
 409                 return( 0 );
 410         }
 411 }
 412 #endif
 413
 414
 415 int includemember( int *mem, int *cand ) // mem in cand
 416 {
 417         int i, j;
 418
 419 #if 0
 420         fprintf( stderr, "mem = " );
 421         for( i=0; mem[i]>-1; i++ )      fprintf( stderr, "%d ", mem[i] );
 422         fprintf( stderr, "\n" );
 423
 424         fprintf( stderr, "cand = " );
 425         for( i=0; cand[i]>-1; i++ )     fprintf( stderr, "%d ", cand[i] );
 426         fprintf( stderr, "\n" );
 427 #endif
 428
 429         for( i=0; mem[i]>-1; i++ )
 430         {
 431                 for( j=0; cand[j]>-1; j++ )
 432                         if( mem[i] == cand[j] ) break;
 433                 if( cand[j] == -1 ) return( 0 );
 434         }
 435 //      fprintf( stderr, "INCLUDED! mem[0]=%d\n", mem[0] );
 436         return( 1 );
 437 }