website/archive/binaries/mac/src/tcoffee/t_coffee_source/util_dp_fasta_nw.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <stdarg.h>
   5 #include <string.h>
   6 #include <ctype.h>
   7
   8 #include "io_lib_header.h"
   9 #include "util_lib_header.h"
  10 #include "define_header.h"
  11 #include "dp_lib_header.h"
  12
  13
  14 int commonsextet( int *table, int *pointt );
  15 void makecompositiontable( int *table, int *pointt );
  16 int *code_seq (char *seq, char *type);
  17 int * makepointtable( int *pointt, int *n, int ktup );
  18
  19 static int tsize;
  20
  21 /**
  22 * calculates the number of common tuples
  23 */
  24 int commonsextet( int *table, int *pointt )
  25 {
  26         int value = 0;
  27         int tmp;
  28         int point;
  29         static int *memo = NULL;
  30         static int *ct = NULL;
  31         static int *cp;
  32
  33         if( !memo )
  34         {
  35                 memo = vcalloc( tsize+1, sizeof( int ) );
  36                 ct = vcalloc( tsize+1, sizeof( int ) );
  37         }
  38
  39         cp = ct;
  40         while( ( point = *pointt++ ) != END_ARRAY )
  41         {
  42           tmp = memo[point]++;
  43           if( tmp < table[point] )
  44             value++;
  45           if( tmp == 0 )
  46             {
  47               *cp++ = point;
  48             }
  49         }
  50         *cp = END_ARRAY;
  51
  52         cp =  ct;
  53         while( *cp != END_ARRAY )
  54                 memo[*cp++] = 0;
  55
  56         return( value );
  57 }
  58
  59 /**
  60 *       calculates how many of each tuple exist
  61 */
  62 void makecompositiontable( int *table, int *pointt )
  63 {
  64         int point;
  65
  66         while( ( point = *pointt++ ) != END_ARRAY )
  67           {
  68             table[point]++;
  69           }
  70 }
  71
  72 int *code_seq (char *seq, char *type)
  73 {
  74   static int *code;
  75   static int *aa, ng;
  76   int a, b, l;
  77
  78
  79   if (!aa)
  80     {
  81       char **gl;
  82       if ( strm (type, "DNA") || strm (type, "RNA"))
  83         {
  84           gl=declare_char (4,5);
  85           sprintf ( gl[ng++], "Aa");
  86           sprintf ( gl[ng++], "Gg");
  87           sprintf ( gl[ng++], "TtUu");
  88           sprintf ( gl[ng++], "Cc");
  89         }
  90       else
  91         {
  92
  93           gl=make_group_aa ( &ng, "mafft");
  94         }
  95       aa=vcalloc ( 256, sizeof (int));
  96       for ( a=0; a<ng; a++)
  97         {
  98           for ( b=0; b< strlen (gl[a]); b++)
  99             {
 100               aa[(int)gl[a][b]]=a;
 101             }
 102         }
 103       free_char (gl, -1);
 104     }
 105
 106
 107   l=strlen (seq);
 108
 109   if ( code) code--;
 110
 111   if ( !code || read_array_size (code, sizeof (int))<(l+2))
 112     {
 113       vfree (code);
 114       code=vcalloc (l+2, sizeof (int));
 115     }
 116   code[0]=ng;
 117   code++;
 118   for (a=0; a<l; a++)
 119     {
 120       code[a]=aa[(int)seq[a]];
 121     }
 122
 123   code[a]=END_ARRAY;
 124   return code;
 125 }
 126
 127
 128 int * makepointtable( int *pointt, int *n, int ktup )
 129 {
 130   int point, a, ng;
 131   register int *p;
 132   static int *prod;
 133
 134   ng=n[-1];
 135
 136   if (!prod)
 137     {
 138       prod=vcalloc ( ktup, sizeof (int));
 139       for ( a=0; a<ktup; a++)
 140         {
 141           prod[ktup-a-1]=(int)pow(n[-1],a);
 142         }
 143     }
 144   p = n;
 145
 146   for (point=0,a=0; a<ktup; a++)
 147     {
 148       point+= *n++ *prod[a];
 149     }
 150
 151   *pointt++ = point;
 152
 153   while( *n != END_ARRAY )
 154     {
 155       point -= *p++ * prod[0];
 156       point *= ng;
 157       point += *n++;
 158       *pointt++ = point;
 159     }
 160   *pointt = END_ARRAY;
 161   return pointt;
 162 }
 163
 164
 165 int ** ktup_dist_mat ( char **seq, int nseq, int ktup, char *type)
 166 {
 167   //Adapted from MAFFT 5: fast ktup
 168   int **pointt,*code=NULL, **pscore;
 169   int i, l, j, minl;
 170   double **mtx, score0;
 171
 172
 173   if (!seq || nseq==0)return NULL;
 174   for (minl=strlen(seq[0]),l=0,i=0;i<nseq; i++)
 175     {
 176       int len;
 177       len=strlen (seq[i]);
 178       minl=MIN(minl, len);
 179       l=MAX(l,len);
 180     }
 181   ktup=MIN(minl, ktup);
 182   pointt=declare_int (nseq, l+1);
 183   mtx=declare_double (nseq, nseq);
 184   pscore=declare_int ( nseq, nseq);
 185
 186   for( i=0; i<nseq; i++ )
 187   {
 188       makepointtable( pointt[i], code=code_seq (seq[i], type),ktup);
 189   }
 190   tsize=(int)pow(code[-1], ktup);
 191
 192   for ( i=0; i<nseq; i++)
 193     {
 194       int *table1;
 195       table1=vcalloc ( tsize,sizeof (int));
 196       makecompositiontable( table1, pointt[i]);
 197       for (j=i; j<nseq; j++)
 198         {
 199           mtx[i][j] = commonsextet( table1, pointt[j] );
 200         }
 201       vfree (table1);
 202     }
 203   for( i=0; i<nseq; i++ )
 204     {
 205       score0 = mtx[i][i];
 206       for( j=0; j<nseq; j++ )
 207         pscore[i][j] = (int)( ( score0 - mtx[MIN(i,j)][MAX(i,j)] ) / score0 * 3 * 10.0 + 0.5 );
 208     }
 209   for( i=0; i<nseq-1; i++ )
 210     for( j=i+1; j<nseq; j++ )
 211       {
 212         pscore[i][j] = pscore[j][i]=100-MIN( pscore[i][j], pscore[j][i] );
 213       }
 214     return pscore;
 215 }
 216
 217
 218 int ** evaluate_diagonals_with_ktup_1 ( Alignment *A, int *ns, int **l_s, Constraint_list *CL,int maximise,int n_groups, char **group_list, int ktup);
 219 int ** evaluate_diagonals_with_ktup_2 ( Alignment *A, int *ns, int **l_s, Constraint_list *CL,int maximise,int n_groups, char **group_list, int ktup);
 220
 221
 222 int ** evaluate_diagonals_for_two_sequences ( char *seq1, char *seq2,int maximise,Constraint_list *CL,int ktup)
 223        {
 224
 225        static int ng;
 226        static char **gl;
 227        static int *ns, **l_s;
 228        Alignment *A;
 229        int **diag;
 230        int in_cl;
 231        char *type;
 232
 233        if (!CL)
 234             {
 235               in_cl=0;
 236
 237               CL=vcalloc ( 1, sizeof (Constraint_list));
 238               CL->maximise=1;
 239               sprintf ( CL->matrix_for_aa_group, "vasiliky");
 240               CL->M=read_matrice ("blosum62mt");
 241               CL->evaluate_residue_pair=evaluate_cdna_matrix_score;
 242               CL->get_dp_cost=slow_get_dp_cost;
 243               type=get_string_type(seq1);
 244
 245               if ( strm (type, "CDNA"))
 246                    CL->evaluate_residue_pair= evaluate_matrix_score;
 247               else if (  strm(type, "PROTEIN"))
 248                    CL->evaluate_residue_pair=evaluate_matrix_score;
 249               else if (  strm (type, "DNA") || strm (type, "RNA"))
 250                    CL->evaluate_residue_pair= evaluate_matrix_score;
 251               vfree(type);
 252             }
 253        else
 254             {
 255               in_cl=1;
 256             }
 257
 258
 259
 260
 261        if ( !gl)
 262          {
 263            gl=make_group_aa (&ng, CL->matrix_for_aa_group);
 264            ns=vcalloc (2, sizeof (int));
 265            ns[0]=ns[1]=1;
 266            l_s=declare_int (2, 2);
 267            l_s[0][0]=0;
 268            l_s[1][0]=1;
 269          }
 270
 271
 272        A=strings2aln (2, "A",seq1,"B", seq2);
 273        ungap(A->seq_al[0]);
 274        ungap(A->seq_al[1]);
 275
 276        CL->S=A->S;
 277
 278        diag=evaluate_diagonals ( A,ns, l_s, CL,maximise, ng, gl, ktup);
 279        free_sequence (A->S, (A->S)->nseq);
 280        free_aln (A);
 281        if (!in_cl)
 282          {
 283           free_int (CL->M, -1);
 284           vfree (CL);
 285          }
 286
 287
 288        return diag;
 289        }
 290
 291
 292 int ** evaluate_diagonals ( Alignment *A, int *ns, int **l_s, Constraint_list *CL,int maximise,int n_groups, char **group_list, int ktup)
 293         {
 294         int **tot_diag;
 295
 296
 297
 298         if      ( CL->residue_index)
 299           {
 300           tot_diag=evaluate_diagonals_with_clist ( A, ns, l_s, CL, maximise,n_groups,group_list, ktup);
 301           }
 302         else if ( CL->use_fragments)
 303             {
 304
 305               tot_diag=evaluate_segments_with_ktup ( A, ns, l_s, CL, maximise,n_groups,group_list, ktup);
 306             }
 307         else
 308           {
 309
 310             tot_diag=evaluate_diagonals_with_ktup ( A, ns, l_s, CL, maximise,n_groups,group_list, ktup);
 311           }
 312
 313         return tot_diag;
 314         }
 315 int ** evaluate_segments_with_ktup ( Alignment *A, int *ns, int **l_s, Constraint_list *CL,int maximise,int n_groups, char **group_list, int ktup)
 316     {
 317    /*
 318     Reads in an alignmnet A, with two groups of sequences marked.
 319     1-Turn each group into a conscensus, using the group list identifier.
 320                -if the group list is left empty original symbols are used
 321     2-hash groupc the two sequences
 322     3-score each diagonal, sort the list and return it (diag_list)
 323    */
 324
 325     char *seq1, *seq2, *alphabet=NULL;
 326     int a,b,l1, l2, n_ktup,pos_ktup1, pos_ktup2, **pos;
 327     int *hasched_seq1, *hasched_seq2,*lu_seq1,*lu_seq2;
 328     int n_diag, **diag, current_diag, **dot_list, n_dots, cost;
 329     int l,delta_diag, delta_res;
 330
 331
 332     pos=aln2pos_simple ( A,-1, ns, l_s);
 333     seq1=aln2cons_seq (A, ns[0], l_s[0], n_groups, group_list);
 334     seq2=aln2cons_seq (A, ns[1], l_s[1], n_groups, group_list);
 335
 336
 337
 338     alphabet=get_alphabet (seq1,alphabet);
 339     alphabet=get_alphabet (seq2,alphabet);
 340
 341
 342
 343     l1=strlen ( seq1);
 344     l2=strlen ( seq2);
 345
 346     n_diag=l1+l2-1;
 347     diag=declare_int ( n_diag+2, 3);
 348     n_ktup=(int)pow ( (double)alphabet[0]+1, (double)ktup);
 349
 350     hasch_seq(seq1, &hasched_seq1, &lu_seq1,ktup, alphabet);
 351     hasch_seq(seq2, &hasched_seq2, &lu_seq2,ktup, alphabet);
 352
 353
 354
 355     /*EVALUATE THE DIAGONALS*/
 356     for ( a=0; a<= n_diag; a++)diag[a][0]=a;
 357     for ( n_dots=0,a=1; a<= n_ktup; a++)
 358         {
 359             pos_ktup1=lu_seq1[a];
 360             while (TRUE)
 361                   {
 362                   if (!pos_ktup1)break;
 363                   pos_ktup2=lu_seq2[a];
 364                   while (pos_ktup2)
 365                             {
 366                             n_dots++;
 367                             pos_ktup2=hasched_seq2[pos_ktup2];
 368                             }
 369                   pos_ktup1=hasched_seq1[pos_ktup1];
 370                   }
 371         }
 372
 373     if ( n_dots==0)
 374        {
 375             vfree (seq1);
 376             vfree (seq2);
 377             vfree (alphabet);
 378             vfree (hasched_seq1);
 379             vfree (hasched_seq2);
 380             vfree (lu_seq1);
 381             vfree (lu_seq2);
 382             free_int (diag, -1);
 383            return evaluate_segments_with_ktup (A,ns,l_s,CL,maximise,n_groups, group_list,ktup-1);
 384        }
 385
 386     dot_list=declare_int ( n_dots,3);
 387
 388     for ( n_dots=0,a=1; a<= n_ktup; a++)
 389         {
 390             pos_ktup1=lu_seq1[a];
 391             while (TRUE)
 392                   {
 393                   if (!pos_ktup1)break;
 394                   pos_ktup2=lu_seq2[a];
 395                   while (pos_ktup2)
 396                             {
 397                             current_diag=(pos_ktup2-pos_ktup1+l1);
 398                             dot_list[n_dots][0]=current_diag;
 399                             dot_list[n_dots][1]=pos_ktup1;
 400                             dot_list[n_dots][2]=pos_ktup2;
 401                             pos_ktup2=hasched_seq2[pos_ktup2];
 402                             n_dots++;
 403                             }
 404                   pos_ktup1=hasched_seq1[pos_ktup1];
 405                   }
 406         }
 407
 408
 409
 410     hsort_list_array ((void **)dot_list, n_dots, sizeof (int), 3, 0, 3);
 411     current_diag= (int)dot_list[0][0];
 412
 413     for ( b=0; b< ktup; b++)diag[current_diag][2]+=(CL->get_dp_cost) ( A, pos, ns[0], l_s[0], dot_list[0][1]+b-1, pos,ns[1], l_s[1], dot_list[0][2]+b-1, CL);
 414
 415
 416     for ( l=0,a=1; a< n_dots; a++)
 417         {
 418
 419             delta_diag=dot_list[a][0]-dot_list[a-1][0];
 420             delta_res =dot_list[a][1]-dot_list[a-1][1];
 421
 422             for ( cost=0, b=0; b< ktup; b++)cost++;
 423
 424             /*=(CL->get_dp_cost) ( A, pos, ns[0], l_s[0], dot_list[a][1]+b-1, pos,ns[1], l_s[1], dot_list[a][2]+b-1, CL);*/
 425
 426
 427
 428             if (delta_diag!=0 || FABS(delta_res)>5)
 429                {
 430
 431                  l=0;
 432                  diag[current_diag][1]=best_of_a_b(diag[current_diag][2], diag[current_diag][1], 1);
 433                  if ( diag[current_diag][2]<0);
 434                  else diag[current_diag][1]= MAX(diag[current_diag][1],diag[current_diag][2]);
 435                  diag[current_diag][2]=0;
 436                  current_diag=dot_list[a][0];
 437                }
 438             l++;
 439             diag[current_diag][2]+=cost;
 440
 441         }
 442     diag[current_diag][1]=best_of_a_b(diag[current_diag][2], diag[current_diag][1], 1);
 443     sort_int (diag+1, 3, 1,0, n_diag-1);
 444
 445
 446     vfree (seq1);
 447     vfree (seq2);
 448     vfree (alphabet);
 449     vfree (hasched_seq1);
 450     vfree (hasched_seq2);
 451     vfree (lu_seq1);
 452     vfree (lu_seq2);
 453     free_int (pos, -1);
 454     free_int (dot_list, -1);
 455     return diag;
 456     }
 457
 458
 459
 460
 461
 462 int ** evaluate_diagonals_with_clist ( Alignment *A, int *ns, int **l_s, Constraint_list *CL,int maximise,int n_groups, char **group_list, int ktup)
 463     {
 464
 465    /*
 466     Reads in an alignmnent A, with two groups of sequences marked.
 467     Weight the diagonals with the values read in the constraint list
 468    */
 469
 470     int l1, l2,n_diag, s1, s2, r1=0, r2=0;
 471     int a, b, c, d;
 472     int **diag;
 473     int **code;
 474     int **pos;
 475     static int *entry;
 476
 477
 478     if ( !entry)entry=vcalloc ( CL->entry_len+1, CL->el_size);
 479     l1=strlen (A->seq_al[l_s[0][0]]);
 480     l2=strlen (A->seq_al[l_s[1][0]]);
 481
 482     n_diag=l1+l2-1;
 483     diag=declare_int ( n_diag+2, 3);
 484     for ( a=0; a<= n_diag; a++)diag[a][0]=a;
 485
 486     A->S=CL->S;
 487     code=seq2aln_pos (A, ns, l_s);
 488     pos =aln2pos_simple ( A,-1, ns, l_s);
 489
 490
 491     for (a=0; a<ns[0]; a++)
 492
 493         {
 494         s1=A->order[l_s[0][a]][0];
 495         for (b=0; b<ns[1]; b++)
 496             {
 497             s2=A->order[l_s[1][b]][0];
 498             for (r1=1; r1<=(A->S)->len[s1]; r1++)
 499               {
 500                 int e;
 501                 for (e=1; e<CL->residue_index[s1][r1][0]; e+=ICHUNK)
 502                   {
 503                     if (CL->residue_index[s1][r1][e+SEQ2]==s2)
 504                       {
 505                         r2=CL->residue_index[s1][r1][e+R2];
 506                         diag[(r2-r1+l1)][1]+=(CL->get_dp_cost) ( A, pos, ns[0], l_s[0],r1-1, pos,ns[1], l_s[1], r2-1, CL);
 507                       }
 508                   }
 509               }
 510             }
 511         }
 512
 513     sort_int (diag+1, 2, 1,0, n_diag-1);
 514
 515     free_int (code,-1);
 516     free_int (pos, -1);
 517     return diag;
 518     }
 519
 520 int * flag_diagonals (int l1, int l2, int **sorted_diag, float T, int window)
 521     {
 522     int a, b, up, low,current_diag,n_diag;
 523     int * slopes;
 524     int *diag_list;
 525     double mean;
 526     double sd;
 527     int use_z_score=1;
 528
 529
 530     n_diag=l1+l2-1;
 531     mean=return_mean_int ( sorted_diag, n_diag+1, 1);
 532
 533     sd  =return_sd_int ( sorted_diag, n_diag+1, 1, (int)mean);
 534
 535     if ( T==0)
 536       {
 537       use_z_score=1;
 538       T=(((double)sorted_diag[n_diag][1]-mean)/sd)/25;
 539       }
 540
 541
 542     diag_list=vcalloc (l1+l2+1, sizeof (int));
 543     slopes=vcalloc ( n_diag+1, sizeof (int));
 544
 545     for ( a=n_diag; a>0; a--)
 546             {
 547             current_diag=sorted_diag[a][0];
 548
 549
 550             if ( !use_z_score && sorted_diag[a][1]>T)
 551                {
 552                    up=MAX(1,current_diag-window);
 553                    low=MIN(n_diag, current_diag+window);
 554                    for ( b=up; b<=low; b++)slopes[b]=1;
 555                }
 556             else if (use_z_score && ((double)sorted_diag[a][1]-mean)/sd>T)
 557               {
 558                 up=MAX(1,current_diag-window);
 559                 low=MIN(n_diag, current_diag+window);
 560                 for ( b=up; b<=low; b++)slopes[b]=1;
 561               }
 562             else break;
 563             }
 564
 565     for ( a=1, b=0; a<=n_diag; a++)
 566         {
 567             b+=slopes[a];
 568         }
 569
 570     slopes[1]=1;
 571     slopes[l1+l2-1]=1;
 572     slopes[l2]=1;
 573     for (a=0; a<= (l1+l2-1); a++)
 574         if ( slopes[a]){diag_list[++diag_list[0]]=a;}
 575
 576     vfree (slopes);
 577
 578     return diag_list;
 579     }
 580 int * extract_N_diag (int l1, int l2, int **sorted_diag, int n_chosen_diag, int window)
 581     {
 582     int a, b, up, low,current_diag,n_diag;
 583     int * slopes;
 584     int *diag_list;
 585
 586
 587     n_diag=l1+l2-1;
 588
 589     diag_list=vcalloc (l1+l2+1, sizeof (int));
 590     slopes=vcalloc ( n_diag+1, sizeof (int));
 591
 592
 593     for ( a=n_diag; a>0 && a>(n_diag-n_chosen_diag); a--)
 594             {
 595             current_diag=sorted_diag[a][0];
 596             up=MAX(1,current_diag-window);
 597             low=MIN(n_diag, current_diag+window);
 598
 599             for ( b=up; b<=low; b++)slopes[b]=1;
 600             }
 601
 602     /*flag bottom right*/
 603     up=MAX(1,1-window);low=MIN(n_diag,1+window);
 604     for ( a=up; a<=low; a++) slopes[a]=1;
 605
 606     /*flag top left */
 607     up=MAX(1,(l1+l2-1)-window);low=MIN(n_diag,(l1+l2-1)+window);
 608     for ( a=up; a<=low; a++) slopes[a]=1;
 609
 610
 611     /*flag MAIN DIAG SEQ1*/
 612     up=MAX(1,l1-window);low=MIN(n_diag,l1+window);
 613     for ( a=up; a<=low; a++) slopes[a]=1;
 614
 615     /*flag MAIN DIAG SEQ2*/
 616     up=MAX(1,l2-window);low=MIN(n_diag,l2+window);
 617     for ( a=up; a<=low; a++) slopes[a]=1;
 618
 619
 620     for (a=0; a<= (l1+l2-1); a++)
 621         if ( slopes[a]){diag_list[++diag_list[0]]=a;}
 622
 623     vfree (slopes);
 624     return diag_list;
 625     }
 626
 627
 628
 629
 630 int cfasta_gotoh_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
 631     {
 632 /*TREATMENT OF THE TERMINAL GAP PENALTIES*/
 633 /*TG_MODE=0---> gop and gep*/
 634 /*TG_MODE=1---> ---     gep*/
 635 /*TG_MODE=2---> ---     ---*/
 636
 637
 638         int maximise;
 639
 640 /*VARIABLES FOR THE MULTIPLE SEQUENCE ALIGNMENT*/
 641         int **tot_diag;
 642
 643         int *diag;
 644         int ktup;
 645         static int n_groups;
 646         static char **group_list;
 647         int score, new_score;
 648         int n_chosen_diag=20;
 649         int step;
 650         int max_n_chosen_diag;
 651         int l1, l2;
 652         /********Prepare Penalties******/
 653
 654
 655         maximise=CL->maximise;
 656         ktup=CL->ktup;
 657
 658         /********************************/
 659
 660
 661
 662
 663         if ( !group_list)
 664            {
 665
 666                group_list=make_group_aa (&n_groups, CL->matrix_for_aa_group);
 667            }
 668
 669         l1=strlen (A->seq_al[l_s[0][0]]);
 670         l2=strlen (A->seq_al[l_s[1][0]]);
 671
 672         if ( !CL->fasta_step)
 673             {
 674             step=MIN(l1,l2);
 675             step=(int) log ((double)MAX(step, 1));
 676             step=MAX(step, 20);
 677             }
 678         else
 679             {
 680                 step=CL->fasta_step;
 681             }
 682
 683
 684         tot_diag=evaluate_diagonals ( A, ns, l_s, CL, maximise,n_groups,group_list, ktup);
 685
 686
 687         max_n_chosen_diag=strlen (A->seq_al[l_s[0][0]])+strlen (A->seq_al[l_s[1][0]])-1;
 688         /*max_n_chosen_diag=(int)log10((double)(l1+l2))*10;*/
 689
 690         n_chosen_diag+=step;
 691         n_chosen_diag=MIN(n_chosen_diag, max_n_chosen_diag);
 692
 693
 694         diag=extract_N_diag (strlen (A->seq_al[l_s[0][0]]),strlen (A->seq_al[l_s[1][0]]), tot_diag, n_chosen_diag, 0);
 695
 696
 697         score    =make_fasta_gotoh_pair_wise ( A, ns, l_s, CL, diag);
 698
 699         new_score=0;
 700         vfree ( diag);
 701
 702
 703         while (new_score!=score && n_chosen_diag< max_n_chosen_diag    )
 704           {
 705
 706
 707             score=new_score;
 708
 709             ungap_sub_aln ( A, ns[0], l_s[0]);
 710             ungap_sub_aln ( A, ns[1], l_s[1]);
 711
 712
 713             n_chosen_diag+=step;
 714             n_chosen_diag=MIN(n_chosen_diag, max_n_chosen_diag);
 715
 716
 717             diag     =extract_N_diag (strlen (A->seq_al[l_s[0][0]]),strlen (A->seq_al[l_s[1][0]]), tot_diag, n_chosen_diag, 0);
 718             new_score=make_fasta_gotoh_pair_wise (  A, ns, l_s, CL, diag);
 719
 720             vfree ( diag);
 721
 722           }
 723
 724         score=new_score;
 725         free_int (tot_diag, -1);
 726
 727         return score;
 728     }
 729
 730 int fasta_gotoh_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
 731     {
 732 /*TREATMENT OF THE TERMINAL GAP PENALTIES*/
 733 /*TG_MODE=0---> gop and gep*/
 734 /*TG_MODE=1---> ---     gep*/
 735 /*TG_MODE=2---> ---     ---*/
 736
 737
 738         int maximise;
 739
 740 /*VARIABLES FOR THE MULTIPLE SEQUENCE ALIGNMENT*/
 741         int **tot_diag;
 742         int *diag;
 743         int ktup;
 744         float diagonal_threshold;
 745         static int n_groups;
 746         static char **group_list;
 747         int score;
 748         /********Prepare Penalties******/
 749
 750
 751         maximise=CL->maximise;
 752         ktup=CL->ktup;
 753         diagonal_threshold=CL->diagonal_threshold;
 754         /********************************/
 755
 756
 757
 758         if ( !group_list)
 759            {
 760                group_list=make_group_aa (&n_groups, CL->matrix_for_aa_group);
 761            }
 762
 763
 764         tot_diag=evaluate_diagonals ( A, ns, l_s, CL, maximise,n_groups,group_list, ktup);
 765
 766         if (  !CL->fasta_step)
 767           {
 768             diag=flag_diagonals (strlen(A->seq_al[l_s[0][0]]),strlen(A->seq_al[l_s[1][0]]), tot_diag,diagonal_threshold,0);
 769           }
 770
 771         else
 772           {
 773
 774             diag=extract_N_diag (strlen (A->seq_al[l_s[0][0]]),strlen (A->seq_al[l_s[1][0]]), tot_diag,CL->fasta_step,0);
 775
 776           }
 777         score=make_fasta_gotoh_pair_wise ( A, ns, l_s, CL, diag);
 778
 779         free_int (tot_diag, -1);
 780         vfree (diag);
 781         return score;
 782     }
 783 int very_fast_gotoh_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
 784     {
 785 /*TREATMENT OF THE TERMINAL GAP PENALTIES*/
 786 /*TG_MODE=0---> gop and gep*/
 787 /*TG_MODE=1---> ---     gep*/
 788 /*TG_MODE=2---> ---     ---*/
 789
 790
 791         int maximise;
 792 /*VARIABLES FOR THE MULTIPLE SEQUENCE ALIGNMENT*/
 793         int **tot_diag;
 794         int *diag;
 795         int ktup;
 796         static int n_groups;
 797         static char **group_list;
 798         int score;
 799         /********Prepare Penalties******/
 800
 801
 802         maximise=CL->maximise;
 803         ktup=CL->ktup;
 804         /********************************/
 805
 806
 807         if ( !group_list)
 808            {
 809
 810                group_list=make_group_aa (&n_groups, CL->matrix_for_aa_group);
 811            }
 812
 813         CL->use_fragments=0;
 814         tot_diag=evaluate_diagonals ( A, ns, l_s, CL, maximise,n_groups,group_list, ktup);
 815
 816         /*Note: 20 diagonals. 5 shadows on each side: tunned on Hom39, 2/2/04 */
 817         diag=extract_N_diag (strlen (A->seq_al[l_s[0][0]]),strlen (A->seq_al[l_s[1][0]]), tot_diag,20,5);
 818         score=make_fasta_gotoh_pair_wise ( A, ns, l_s, CL, diag);
 819         free_int (tot_diag, -1);
 820         vfree (diag);
 821         return score;
 822     }
 823 int make_fasta_gotoh_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL, int *diag)
 824     {
 825 /*TREATMENT OF THE TERMINAL GAP PENALTIES*/
 826 /*TG_MODE=0---> gop and gep*/
 827 /*TG_MODE=1---> ---     gep*/
 828       /*TG_MODE=2---> ---     ---*/
 829
 830
 831         int TG_MODE, gop, l_gop, gep,l_gep, maximise;
 832
 833 /*VARIABLES FOR THE MULTIPLE SEQUENCE ALIGNMENT*/
 834         int a, b,c,k, t;
 835         int l1, l2,eg, ch, sub,score=0, last_i=0, last_j=0, i, delta_i, j, pos_j, ala, alb, LEN, n_diag, match1, match2;
 836         int su, in, de, tr;
 837
 838         int **C, **D, **I, **trace, **pos0, **LD;
 839         int lenal[2], len;
 840         char *buffer, *char_buf;
 841         char **aln, **al;
 842
 843         /********Prepare Penalties******/
 844         gop=CL->gop*SCORE_K;
 845         gep=CL->gep*SCORE_K;
 846         TG_MODE=CL->TG_MODE;
 847         maximise=CL->maximise;
 848
 849
 850         /********************************/
 851
 852
 853         n_diag=diag[0];
 854
 855
 856
 857        l1=lenal[0]=strlen (A->seq_al[l_s[0][0]]);
 858        l2=lenal[1]=strlen (A->seq_al[l_s[1][0]]);
 859
 860        if ( getenv ("DEBUG_TCOFFEE"))fprintf ( stderr, "\n\tNdiag=%d%%  ", (diag[0]*100)/(l1+l2));
 861
 862         /*diag:
 863           diag[1..n_diag]--> flaged diagonal in order;
 864           diag[0]=0--> first diagonal;
 865           diag[n_diag+1]=l1+l2-1;
 866         */
 867
 868         /*numeration of the diagonals strats from the bottom right [1...l1+l2-1]*/
 869         /*sequence s1 is vertical and seq s2 is horizontal*/
 870         /*D contains the best Deletion  in S2==>comes from diagonal N+1*/
 871         /*I contains the best insertion in S2=> comes from diagonal N-1*/
 872
 873
 874
 875
 876
 877        C=declare_int (lenal[0]+lenal[1]+1, n_diag+2);
 878        D=declare_int (lenal[0]+lenal[1]+1, n_diag+2);
 879        LD=declare_int (lenal[0]+lenal[1]+1, n_diag+2);
 880        I=declare_int (lenal[0]+lenal[1]+1, n_diag+2);
 881        trace=declare_int (lenal[0]+lenal[1]+1, n_diag+2);
 882
 883
 884        al=declare_char (2,lenal[0]+lenal[1]+lenal[1]+1);
 885
 886        len= MAX(lenal[0],lenal[1])+1;
 887        buffer=vcalloc ( 2*len, sizeof (char));
 888        char_buf= vcalloc (2*len, sizeof (char));
 889
 890        pos0=aln2pos_simple ( A,-1, ns, l_s);
 891        C[0][0]=0;
 892
 893        t=(TG_MODE==0)?gop:0;
 894        for ( j=1; j<= n_diag; j++)
 895             {
 896                 l_gop=(TG_MODE==0)?gop:0;
 897                 l_gep=(TG_MODE==2)?0:gep;
 898
 899
 900
 901                 if ( (diag[j]-lenal[0])<0 )
 902                     {
 903                     trace[0][j]=UNDEFINED;
 904                     continue;
 905                     }
 906                 C[0][j]=(diag[j]-lenal[0])*l_gep +l_gop;
 907                 D[0][j]=(diag[j]-lenal[0])*l_gep +l_gop+gop;
 908             }
 909        D[0][j]=D[0][j-1]+gep;
 910
 911
 912        t=(TG_MODE==0)?gop:0;
 913        for ( i=1; i<=lenal[0]; i++)
 914            {
 915                 l_gop=(TG_MODE==0)?gop:0;
 916                 l_gep=(TG_MODE==2)?0:gep;
 917
 918                 C[i][0]=C[i][n_diag+1]=t=t+l_gep;
 919                 I[i][0]=D[i][n_diag+1]=t+    gop;
 920
 921                 for ( j=1; j<=n_diag; j++)
 922                     {
 923                         C[i][j]=C[i][0];
 924                         D[i][j]=I[i][j]=I[i][0];
 925                     }
 926
 927                 for (eg=0, j=1; j<=n_diag; j++)
 928                     {
 929
 930                         pos_j=diag[j]-lenal[0]+i;
 931                         if (pos_j<=0 || pos_j>l2 )
 932                             {
 933                             trace[i][j]=UNDEFINED;
 934                             continue;
 935                             }
 936                         sub=(CL->get_dp_cost) ( A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],pos_j-1, CL );
 937
 938                     /*1 identify the best insertion in S2:*/
 939                         l_gop=(i==lenal[0])?((TG_MODE==0)?gop:0):gop;
 940                         l_gep=(i==lenal[0])?((TG_MODE==2)?0:gep):gep;
 941                         len=(j==1)?0:(diag[j]-diag[j-1]);
 942                         if ( a_better_than_b(I[i][j-1], C[i][j-1]+l_gop, maximise))eg++;
 943                         else eg=1;
 944                         I[i][j]=best_of_a_b (I[i][j-1], C[i][j-1]+l_gop, maximise)+len*l_gep;
 945
 946                     /*2 Identify the best deletion in S2*/
 947                         l_gop=(pos_j==lenal[1])?((TG_MODE==0)?gop:0):gop;
 948                         l_gep=(pos_j==lenal[1])?((TG_MODE==2)?0:gep):gep;
 949
 950                         len=(j==n_diag)?0:(diag[j+1]-diag[j]);
 951                         delta_i=((i-len)>0)?(i-len):0;
 952
 953                         if ( a_better_than_b(D[delta_i][j+1],C[delta_i][j+1]+l_gop, maximise)){LD[i][j]=LD[delta_i][j+1]+1;}
 954                         else {LD[i][j]=1;}
 955                         D[i][j]=best_of_a_b (D[delta_i][j+1],C[delta_i][j+1]+l_gop, maximise)+len*l_gep;
 956
 957
 958                         /*Identify the best way*/
 959                         /*
 960                         score=C[i][j]=best_int ( 3, maximise, &fop, I[i][j], C[i-1][j]+sub, D[i][j]);
 961                         fop-=1;
 962                         if ( fop<0)trace[i][j]=fop*eg;
 963                         else if ( fop>0 ) {trace[i][j]=fop*LD[i][j];}
 964                         else if ( fop==0) trace[i][j]=0;
 965                         */
 966
 967                         su=C[i-1][j]+sub;
 968                         in=I[i][j];
 969                         de=D[i][j];
 970
 971                         /*HERE ("%d %d %d", su, in, de);*/
 972                         if (su>=in && su>=de)
 973                           {
 974                             score=su;
 975                             tr=0;
 976                           }
 977                         else if (in>=de)
 978                           {
 979                             score=in;
 980                             tr=-eg;
 981                           }
 982                         else
 983                           {
 984                             score=de;
 985                             tr=LD[i][j];
 986                           }
 987                         trace[i][j]=tr;
 988                         C[i][j]=score;
 989
 990
 991                         last_i=i;
 992                         last_j=j;
 993                     }
 994             }
 995
 996
 997        /*
 998                     [0][Positive]
 999                      ^     ^
1000                      |    /
1001                      |   /
1002                      |  /
1003                      | /
1004                      |/
1005        [Neg]<-------[*]
1006         */
1007
1008
1009         i=last_i;
1010         j=last_j;
1011
1012
1013
1014         ala=alb=0;
1015         match1=match2=0;
1016         while (!(match1==l1 && match2==l2))
1017               {
1018
1019
1020                   if ( match1==l1)
1021                      {
1022                          len=l2-match2;
1023                          for ( a=0; a< len; a++)
1024                              {
1025                              al[0][ala++]=0;
1026                              al[1][alb++]=1;
1027                              match2++;
1028                              }
1029                          k=0;
1030                          break;
1031
1032                          /*k=-(j-1);*/
1033
1034                      }
1035                   else if ( match2==l2)
1036                      {
1037                          len=l1-match1;
1038                          for ( a=0; a< len; a++)
1039                              {
1040                              al[0][ala++]=1;
1041                              al[1][alb++]=0;
1042                              match1++;
1043                              }
1044                          k=0;
1045                          break;
1046                          /*k= n_diag-j;*/
1047                      }
1048                   else
1049                       {
1050                           k=trace[i][j];
1051                       }
1052
1053
1054                   if ( k==0)
1055                              {
1056                                  if ( match2==l2 || match1==l1);
1057                                  else
1058                                     {
1059
1060                                     al[0][ala++]=1;
1061                                     al[1][alb++]=1;
1062                                     i--;
1063                                     match1++;
1064                                     match2++;
1065                                     }
1066                              }
1067                   else if ( k>0)
1068                              {
1069
1070                              len=diag[j+k]-diag[j];
1071                              for ( a=0; a<len; a++)
1072                                  {
1073                                      if ( match1==l1)break;
1074                                      al[0][ala++]=1;
1075                                      al[1][alb++]=0;
1076                                      match1++;
1077                                  }
1078                              i-=len;
1079                              j+=k;
1080                              }
1081                   else if ( k<0)
1082                              {
1083                              k*=-1;
1084                              len=diag[j]-diag[j-k];
1085                              for ( a=0; a<len; a++)
1086                                  {
1087                                      if ( match2==l2)break;
1088                                      al[0][ala++]=0;
1089                                      al[1][alb++]=1;
1090                                      match2++;
1091                                  }
1092
1093
1094                              j-=k;
1095                              }
1096               }
1097
1098         LEN=ala;
1099         c=LEN-1;
1100         invert_list_char ( al[0], LEN);
1101         invert_list_char ( al[1], LEN);
1102         if ( A->declared_len<=LEN)A=realloc_aln2  ( A,A->max_n_seq, 2*LEN);
1103         aln=A->seq_al;
1104
1105         for ( c=0; c< 2; c++)
1106             {
1107             for ( a=0; a< ns[c]; a++)
1108                 {
1109                 ch=0;
1110                 for ( b=0; b< LEN; b++)
1111                     {
1112                     if (al[c][b]==1)
1113                         char_buf[b]=aln[l_s[c][a]][ch++];
1114                     else
1115                         char_buf[b]='-';
1116                    }
1117                 char_buf[b]='\0';
1118                 sprintf (aln[l_s[c][a]],"%s", char_buf);
1119                 }
1120              }
1121
1122
1123         A->len_aln=LEN;
1124         A->nseq=ns[0]+ns[1];
1125
1126         free_int (pos0, -1);
1127         free_int (C, -1);
1128         free_int (D, -1);
1129         free_int (I, -1);
1130         free_int (trace, -1);
1131         free_int (LD, -1);
1132         free_char ( al, -1);
1133         vfree(buffer);
1134         vfree(char_buf);
1135
1136
1137         return score;
1138     }
1139
1140 int hasch_seq(char *seq, int **hs, int **lu,int ktup,char *alp)
1141     {
1142         static int a[10];
1143
1144         int i,j,l,limit,code,flag;
1145         char residue;
1146
1147         int alp_lu[10000];
1148         int alp_size;
1149
1150         alp_size=alp[0];
1151         alp++;
1152
1153
1154
1155         for ( i=0; i< alp_size; i++)
1156             {
1157               alp_lu[(int)alp[i]]=i;
1158             }
1159
1160
1161
1162         l=strlen (seq);
1163         limit = (int)   pow((double)(alp_size+1),(double)ktup);
1164         hs[0]=vcalloc ( l+1,sizeof (int));
1165         lu[0]=vcalloc ( limit+1, sizeof(int));
1166
1167
1168         if ( l==0)myexit(EXIT_FAILURE);
1169
1170         for (i=1;i<=ktup;i++)
1171            a[i] = (int) pow((double)(alp_size+1),(double)(i-1));
1172
1173
1174         for(i=1;i<=(l-ktup+1);++i)
1175                 {
1176                 code=0;
1177                 flag=FALSE;
1178                 for(j=1;j<=ktup;++j)
1179                    {
1180                    if (is_gap(seq[i+j-2])){flag=TRUE;break;}
1181                    else residue=alp_lu[(int)seq[i+j-2]];
1182                    code+=residue*a[j];
1183                    }
1184
1185                 if ( flag)continue;
1186                 ++code;
1187
1188                 if (lu[0][code])hs[0][i]=lu[0][code];
1189                 lu[0][code]=i;
1190                 }
1191         return 0;
1192     }
1193
1194
1195
1196 /*********************************************************************/
1197 /*                                                                   */
1198 /*                         KTUP_DP                                   */
1199 /*                                                                   */
1200 /*                                                                   */
1201 /*********************************************************************/
1202
1203 /**************Hasch DAta Handling*******************************************************/
1204
1205 struct Hasch_data * free_ktup_hasch_data (struct Hasch_data *d);
1206 struct Hasch_data * declare_ktup_hasch_data (struct Hasch_entry *e);
1207 struct Hasch_data * allocate_ktup_hasch_data (struct Hasch_data *e, int action);
1208
1209 struct Hasch_data
1210 {
1211  int *list;
1212 };
1213 typedef struct Hasch_data Hasch_data;
1214 struct Hasch_data * free_ktup_hasch_data (struct Hasch_data *d)
1215 {
1216   return allocate_ktup_hasch_data (d, FREE);
1217 }
1218 struct Hasch_data * declare_ktup_hasch_data (struct Hasch_entry *e)
1219 {
1220   e->data=allocate_ktup_hasch_data (NULL,DECLARE);
1221   return e->data;
1222 }
1223
1224 struct Hasch_data * allocate_ktup_hasch_data (struct Hasch_data *e, int action)
1225 {
1226   static struct Hasch_data **heap;
1227   static int heap_size, free_heap, a;
1228
1229   if ( action == 100)
1230     {
1231       fprintf ( stderr, "\nHeap size: %d, Free Heap: %d", heap_size, free_heap);
1232       return NULL;
1233     }
1234   else if ( action==DECLARE)
1235     {
1236       if ( free_heap==0)
1237         {
1238           free_heap=100;
1239           heap_size+=free_heap;
1240           heap=vrealloc (heap,heap_size*sizeof (struct Hasch_entry *));
1241           for ( a=0; a<free_heap; a++)
1242             {
1243               (heap[a])=vcalloc ( 1, sizeof ( struct Hasch_entry *));
1244               (heap[a])->list=vcalloc ( 10, sizeof (int));
1245               (heap[a])->list[0]=10;
1246             }
1247         }
1248       return heap[--free_heap];
1249     }
1250   else if ( action==FREE)
1251     {
1252       heap[free_heap++]=e;
1253       e->list[1]=0;
1254       return NULL;
1255     }
1256   return NULL;
1257 }
1258
1259
1260 /**************Hasch DAta Handling*******************************************************/
1261
1262 int precomputed_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
1263     {
1264       int l1, l2, a, b, c;
1265       int nid=0, npos=0, id;
1266       int r1, r2, s1, s2;
1267
1268       l1=strlen(A->seq_al[l_s[0][0]]);
1269       l2=strlen(A->seq_al[l_s[1][0]]);
1270       if (l1!=l2)
1271         {
1272           fprintf ( stderr, "\nERROR: improper use of the function precomputed pairwise:[FATAL:%s]", PROGRAM);
1273           crash ("");
1274         }
1275       else if ( l1==0)
1276         {
1277           A->score_aln=A->score=0;
1278           return 0;
1279         }
1280
1281       for (npos=0, nid=0, a=0; a< ns[0]; a++)
1282         {
1283           s1=l_s[0][a];
1284
1285           for (b=0; b< ns[1]; b++)
1286             {
1287               s2=l_s[1][b];
1288               for ( c=0; c<l1; c++)
1289                 {
1290                 r1=A->seq_al[s1][c];
1291                 r2=A->seq_al[s2][c];
1292                 if ( is_gap(r1) || is_gap(r2));
1293                 else
1294                   {
1295                     npos++;
1296                     nid+=(r1==r2);
1297                   }
1298                 }
1299             }
1300         }
1301       id=(npos==0)?0:((nid*100)/npos);
1302       A->score=A->score_aln=id;
1303       return A->score;
1304     }
1305 int ktup_comparison_str ( char *seq1, char *seq2, const int ktup);
1306 int ktup_comparison_hasch ( char *i_seq1, char *i_seq2, const int ktup);
1307 int ktup_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
1308     {
1309       static char **gl;
1310       static int ng;
1311       char *seq1;
1312       char *seq2;
1313
1314       int min_len=10;
1315
1316
1317
1318       if ( !gl)
1319         gl=make_group_aa (&ng, "vasiliky");
1320
1321
1322       if ( ns[0]>1)seq1=sub_aln2cons_seq_mat (A, ns[0], l_s[0],"blosum62mt");
1323       else
1324         {
1325           seq1=vcalloc ( strlen (A->seq_al[l_s[0][0]])+1, sizeof (char));
1326           sprintf ( seq1, "%s",A->seq_al[l_s[0][0]]);
1327         }
1328       if ( ns[1]>1)seq2=sub_aln2cons_seq_mat (A, ns[1], l_s[1],"blosum62mt");
1329       else
1330         {
1331           seq2=vcalloc ( strlen (A->seq_al[l_s[1][0]])+1, sizeof (char));
1332           sprintf ( seq2, "%s",A->seq_al[l_s[1][0]]);
1333         }
1334
1335       if ( strlen (seq1)<min_len || strlen (seq2)<min_len)
1336         {
1337           Alignment *B;
1338
1339           ungap(seq1); ungap(seq2);
1340           B=align_two_sequences ( seq1, seq2, "blosum62mt",-10, -1, "myers_miller_pair_wise");
1341           A->score=A->score_aln=aln2sim(B, "idmat");
1342           free_aln (B);
1343           return A->score;
1344         }
1345       else
1346         {
1347
1348           string_convert (seq1, ng, gl);
1349           string_convert (seq2, ng, gl);
1350           A->score=A->score_aln=ktup_comparison (seq1,seq2, CL->ktup);
1351         }
1352
1353       vfree (seq1); vfree (seq2);
1354       return A->score;
1355     }
1356 int ktup_comparison( char *seq2, char *seq1, const int ktup)
1357 {
1358   return ktup_comparison_hasch ( seq2, seq1, ktup);
1359 }
1360 int ktup_comparison_str ( char *seq2, char *seq1, const int ktup)
1361 {
1362   int a,l1, l2,c1, c2, end, start;
1363   char *s1, *s2;
1364   double score=0;
1365   int max_dist=-1;
1366
1367   if ( max_dist==-1)max_dist=MAX((strlen (seq1)),(strlen (seq2)));
1368   l1=strlen (seq1)-ktup;
1369   l2=strlen (seq2);
1370
1371
1372   for ( a=0; a< l1; a++)
1373     {
1374       c1=seq1[a+ktup];seq1[a+ktup]='\0';
1375       s1=seq1+a;
1376
1377       start=((a-max_dist)<0)?0:a-max_dist;
1378       end=((a+max_dist)>=l2)?l2:a+max_dist;
1379
1380       c2=seq2[end];seq2[end]='\0';
1381       s2=seq2+start;
1382
1383       score+=(strstr(s2, s1)!=NULL)?1:0;
1384
1385       seq1[a+ktup]=c1;
1386       seq2[end]=c2;
1387     }
1388   score/=(l1==0)?1:l1;
1389   score=((log(0.1+score)-log(0.1))/(log(1.1)-log(0.1)));
1390
1391   return score*100;
1392
1393 }
1394 int ktup_comparison_hasch ( char *i_seq1, char *i_seq2, const int ktup)
1395 {
1396   /*Ktup comparison adapted from Rob Edgar, NAR, vol32, No1, 381, 2004*/
1397   /*1: hasch sequence 1
1398     2: Count the number of seq2 ktup found in seq1
1399   */
1400
1401   char c;
1402   int key;
1403
1404   static HaschT*H1;
1405   static char *pseq;
1406   Hasch_entry *e;
1407   char *s;
1408   int l, ls;
1409   int p, a, max_dist=-1;
1410   double score=0;
1411
1412
1413
1414   if (!strm (i_seq1, pseq))
1415     {
1416       if (H1)
1417         {
1418           hdestroy (H1, declare_ktup_hasch_data, free_ktup_hasch_data);
1419           string2key (NULL, NULL);
1420         }
1421       H1=hasch_sequence ( i_seq1, ktup);
1422       vfree (pseq);pseq=vcalloc ( strlen (i_seq1)+1, sizeof (char));
1423       sprintf ( pseq, "%s", i_seq1);
1424     }
1425
1426   ls=l=strlen (i_seq2);
1427   s=i_seq2;
1428   p=0;
1429   while (ls>ktup)
1430     {
1431       c=s[ktup];s[ktup]='\0';
1432       key=string2key (s, NULL);
1433       e=hsearch (H1,key,FIND, declare_ktup_hasch_data, free_ktup_hasch_data);
1434
1435       if ( e==NULL);
1436       else if ( max_dist==-1)score++;
1437       else
1438         {
1439           for ( a=1; a<=(e->data)->list[1]; a++)
1440             if (FABS((p-(e->data)->list[a]))<=max_dist)
1441               {score++; break;}
1442         }
1443       s[ktup]=c;s++;p++;ls--;
1444     }
1445   score/=(l-ktup);
1446   score=(log(0.1+score)-log(0.1))/(log(1.1)-log(0.1));
1447
1448   if ( score>100) score=100;
1449   return (int)(score*100);
1450 }
1451
1452 HaschT* hasch_sequence ( char *seq1, int ktup)
1453 {
1454   char c;
1455   int key, offset=0, ls;
1456   HaschT *H;
1457   Hasch_entry *e;
1458
1459   H=hcreate ( strlen (seq1), declare_ktup_hasch_data, free_ktup_hasch_data);
1460   ls=strlen (seq1);
1461   while (ls>=(ktup))
1462     {
1463       c=seq1[ktup];seq1[ktup]='\0';
1464       key=string2key (seq1, NULL);
1465       e=hsearch (H,key,FIND, declare_ktup_hasch_data, free_ktup_hasch_data);
1466
1467       if (e==NULL)
1468         {
1469          e=hsearch (H,key,ADD,declare_ktup_hasch_data,free_ktup_hasch_data);
1470          (e->data)->list[++(e->data)->list[1]+1]=offset;
1471         }
1472       else
1473         {
1474           if ((e->data)->list[0]==((e->data)->list[1]+2)){(e->data)->list[0]+=10;(e->data)->list=vrealloc ((e->data)->list,(e->data)->list[0]*sizeof (int));}
1475           (e->data)->list[++(e->data)->list[1]+1]=offset;
1476         }
1477        seq1[ktup]=c;seq1++;ls--;
1478        offset++;
1479     }
1480   return H;
1481 }
1482
1483
1484
1485 char *dayhoff_translate (char *seq1)
1486 {
1487 int l, a, c;
1488 l=strlen (seq1);
1489  for ( a=0; a< l; a++)
1490   {
1491     c=tolower(seq1[a]);
1492     if ( strchr ("agpst", c))seq1[a]='a';
1493     else if (strchr ("denq", c))seq1[a]='d';
1494     else if (strchr ("fwy", c))seq1[a]='f';
1495     else if (strchr ("hkr", c))seq1[a]='h';
1496     else if (strchr ("ilmv", c))seq1[a]='i';
1497   }
1498 return seq1;
1499 }
1500
1501 int ** evaluate_diagonals_with_ktup ( Alignment *A, int *ns, int **l_s, Constraint_list *CL,int maximise,int n_groups, char **group_list, int ktup)
1502 {
1503   /*Ktup comparison as in Rob Edgar, NAR, vol32, No1, 381, 2004*/
1504   char character;
1505   int key,ls;
1506   HaschT*H1, *H2;
1507   Hasch_entry *e1, *e2;
1508   char *s, *sb, *seq1, *seq2;
1509   int l1, l2;
1510   int score=0;
1511   int **diag,n_diag, ktup1, ktup2,a,b,c,d, **pos;
1512   int n_dots=0;
1513
1514   pos=aln2pos_simple ( A,-1, ns, l_s);
1515
1516   seq1=aln2cons_maj (A, ns[0], l_s[0], n_groups, group_list);
1517   seq2=aln2cons_maj (A, ns[1], l_s[1], n_groups, group_list);
1518   l1=strlen (seq1);
1519   l2=strlen (seq2);
1520   n_diag=l1+l2-1;
1521
1522
1523   diag=declare_int (n_diag+2, 3);
1524   for ( a=0; a<n_diag+2; a++)diag[a][0]=a;
1525
1526   H1=hasch_sequence ( seq1, ktup);
1527   H2=hasch_sequence ( seq2, ktup);
1528   s=sb=vcalloc (strlen (seq1)+strlen (seq2)+1, sizeof (char));
1529   sprintf (s, "%s%s", seq1, seq2);
1530
1531   ls=strlen(s);
1532   while (ls>=(ktup))
1533     {
1534       character=s[ktup];s[ktup]='\0';
1535       key=string2key (s, NULL);
1536       e1=hsearch (H1,key,FIND,declare_ktup_hasch_data, free_ktup_hasch_data);
1537       e2=hsearch (H2,key,FIND,declare_ktup_hasch_data, free_ktup_hasch_data);
1538       if ( !e2 || !e1);
1539       else
1540         {
1541
1542           for (b=2; b<(e1->data)->list[1]+2; b++)
1543             for (c=2; c<(e2->data)->list[1]+2; c++)
1544               {
1545
1546                 ktup1=(e1->data)->list[b];
1547                 ktup2=(e2->data)->list[c];
1548                 diag[(ktup2-ktup1)+l1][2]++;
1549                 for (score=0, d=0; d<ktup; d++)
1550                   score+=(CL->get_dp_cost) ( A, pos, ns[0], l_s[0], ktup1+d, pos,ns[1], l_s[1], ktup2+d, CL);
1551                 diag[(ktup2-ktup1)+l1][1]+=score;
1552                 n_dots++;
1553               }
1554           (e1->data)->list[1]=(e2->data)->list[1]=0;
1555         }
1556       s[ktup]=character;s++;ls--;
1557     }
1558
1559   sort_int (diag+1, 2, 1,0,n_diag-1);
1560
1561   hdestroy (H1,declare_ktup_hasch_data, free_ktup_hasch_data); hdestroy (H2,declare_ktup_hasch_data, free_ktup_hasch_data);
1562   vfree (seq1); vfree (seq2);vfree (sb);free_int (pos, -1);
1563   return diag;
1564 }
1565  /*********************************************************************/
1566 /*                                                                   */
1567 /*                         OLD FUNCTIONS                              */
1568 /*                                                                   */
1569 /*                                                                   */
1570 /*********************************************************************/
1571 int ** evaluate_diagonals_with_ktup_1 ( Alignment *A, int *ns, int **l_s, Constraint_list *CL,int maximise,int n_groups, char **group_list, int ktup)
1572     {
1573    /*
1574     Reads in an alignmnent A, with two groups of sequences marked.
1575     1-Turn each group into a conscensus, using the group list identifier.
1576                -if the group list is left empty original symbols are used
1577     2-hasch the two sequences
1578     3-score each diagonal, sort the list and return it (diag_list)
1579
1580         diag_list:
1581
1582    */
1583
1584     char *seq1, *seq2, *alphabet=NULL;
1585     int a,b,l1, l2, n_ktup,pos_ktup1, pos_ktup2, **pos;
1586     int *hasched_seq1, *hasched_seq2,*lu_seq1,*lu_seq2;
1587     int n_diag, **diag, current_diag, n_dots;
1588     static char *buf;
1589     pos=aln2pos_simple ( A,-1, ns, l_s);
1590
1591
1592     seq1=aln2cons_seq (A, ns[0], l_s[0], n_groups, group_list);
1593     seq2=aln2cons_seq (A, ns[1], l_s[1], n_groups, group_list);
1594
1595
1596
1597
1598     alphabet=get_alphabet (seq1,alphabet);
1599     alphabet=get_alphabet (seq2,alphabet);
1600
1601     l1=strlen ( seq1);
1602     l2=strlen ( seq2);
1603
1604     n_diag=l1+l2-1;
1605     diag=declare_int ( n_diag+2, 3);
1606     n_ktup=(int)pow ( (double)alphabet[0]+1, (double)ktup);
1607
1608
1609     hasch_seq(seq1, &hasched_seq1, &lu_seq1,ktup, alphabet);
1610     hasch_seq(seq2, &hasched_seq2, &lu_seq2,ktup, alphabet);
1611
1612
1613
1614
1615     /*EVALUATE THE DIAGONALS*/
1616     for ( a=0; a<= n_diag; a++)diag[a][0]=a;
1617     for ( n_dots=0,a=1; a<= n_ktup; a++)
1618         {
1619             pos_ktup1=lu_seq1[a];
1620             while (TRUE)
1621                   {
1622                   if (!pos_ktup1)break;
1623                   pos_ktup2=lu_seq2[a];
1624                   while (pos_ktup2)
1625                             {
1626                             current_diag=(pos_ktup2-pos_ktup1+l1);
1627                             for ( b=0; b< ktup; b++)
1628                                 {
1629                                     diag[current_diag][1]+=(CL->get_dp_cost) ( A, pos, ns[0], l_s[0], pos_ktup1+b-1, pos,ns[1], l_s[1], pos_ktup2+b-1, CL);
1630                                     n_dots++;
1631
1632                                 }
1633                             diag[current_diag][2]++;
1634                             pos_ktup2=hasched_seq2[pos_ktup2];
1635                             }
1636                   pos_ktup1=hasched_seq1[pos_ktup1];
1637                   }
1638
1639         }
1640     if ( n_dots==0)
1641        {
1642            if ( !buf)
1643                {
1644                buf=vcalloc ( 30, sizeof (30));
1645                sprintf ( buf, "abcdefghijklmnopqrstuvwxyz");
1646                }
1647             vfree ( hasched_seq1);
1648             vfree ( hasched_seq2);
1649             vfree (lu_seq1);
1650             vfree (lu_seq2);
1651            return evaluate_diagonals_with_ktup ( A,ns,l_s, CL,maximise,1,&buf,1);
1652        }
1653
1654
1655     sort_int (diag+1, 2, 1,0, n_diag-1);
1656     vfree (seq1);
1657     vfree (seq2);
1658     vfree (alphabet);
1659     vfree ( hasched_seq1);
1660     vfree ( hasched_seq2);
1661     vfree (lu_seq1);
1662     vfree (lu_seq2);
1663     free_int (pos, -1);
1664     return diag;
1665     }
1666 /////////////////////////////////////////////////////////////////
1667
1668 Constraint_list * hasch2constraint_list (Sequence*S, Constraint_list *CL)
1669 {
1670   int a,b,c, n;
1671   SeqHasch h,*H=NULL;
1672   int *entry;
1673   int ktup=2;
1674
1675
1676   entry=vcalloc ( CL->entry_len+1, sizeof (int));
1677
1678   for (a=0; a<S->nseq; a++)
1679     {
1680       H=seq2hasch (a, S->seq[a],ktup,H);
1681     }
1682
1683   n=1;
1684   while (H[n])
1685     {
1686       h=H[n];
1687
1688       for (a=0; a<h->n-2; a+=2)
1689         {
1690           for (b=a+2; b<h->n; b+=2)
1691             {
1692
1693               if (h->l[a]==h->l[b])continue;
1694               else
1695                 {
1696                   for (c=0; c<ktup; c++)
1697                     {
1698                       entry[SEQ1]=h->l[a];
1699                       entry[SEQ2]=h->l[b];
1700                       entry[R1]=h->l[a+1]+c;
1701                       entry[R2]=h->l[b+1]+c;
1702                       entry[WE]=100;
1703                       add_entry2list (entry,CL);
1704                     }
1705                 }
1706             }
1707         }
1708       n++;
1709     }
1710
1711   return CL;
1712 }
1713 SeqHasch *cleanhasch       (SeqHasch *H)
1714 {
1715   int n=1;
1716   SeqHasch *N;
1717   N=vcalloc (2, sizeof (SeqHasch));
1718   N[0]=H[0];
1719
1720   while (H[n])
1721     {
1722       (H[n])->n=0;
1723       vfree ((H[n])->l);
1724       (H[n])->l=NULL;
1725       n++;
1726     }
1727   vfree (H);
1728   return N;
1729 }
1730 int hasch2sim        (SeqHasch *H, int nseq)
1731 {
1732   int n=1;
1733
1734   int a,cs, ps, ns;
1735   int id=0, tot=0;
1736
1737   while (H[n])
1738     {
1739       for (ps=-1,ns=0,a=0; a<(H[n])->n; a+=2)
1740         {
1741           //HERE ("%d--[%d %d]",n, (H[n])->l[a], (H[n])->l[a+1]);
1742           cs=(H[n])->l[a];
1743           if (cs!=ps)ns++;
1744           ps=cs;
1745         }
1746       n++;
1747       if (ns==nseq)id++;
1748       tot++;
1749     }
1750
1751   return (id*MAXID)/tot;
1752 }
1753 SeqHasch * seq2hasch (int i,char *seq, int ktup, SeqHasch *H)
1754 {
1755   int a,b,l, n=0;
1756   SeqHasch h;
1757
1758
1759   if (!H)
1760     {
1761       H=vcalloc (2, sizeof (SeqHasch));
1762       H[0]=vcalloc (1, sizeof (hseq));
1763       n=1;
1764     }
1765   else
1766     {
1767       n=0;
1768       while (H[++n]);
1769     }
1770
1771   l=strlen (seq);
1772   for (a=0; a<l-ktup; a++)
1773     {
1774       h=H[0];
1775       for (b=a; b<a+ktup; b++)
1776         {
1777           char r;
1778           r=seq[b];
1779           if (!h->hl[r])  h->hl[r]=vcalloc (1, sizeof (hseq));
1780           h=h->hl[r];
1781         }
1782       if (!h->l)
1783         {
1784
1785           h->n=2;
1786           h->l=vcalloc (2, sizeof (int));
1787           H=vrealloc (H,(n+2)*sizeof (SeqHasch));
1788           H[n]=h;
1789           n++;
1790         }
1791       else
1792         {
1793           h->n+=2;
1794           h->l=vrealloc (h->l, (h->n)*sizeof (int));
1795         }
1796
1797       h->l[h->n-2]=i;
1798       h->l[h->n-1]=a;
1799     }
1800   return H;
1801 }
1802
1803 /******************************COPYRIGHT NOTICE*******************************/
1804 /*© Centro de Regulacio Genomica */
1805 /*and */
1806 /*Cedric Notredame */
1807 /*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
1808 /*All rights reserved.*/
1809 /*This file is part of T-COFFEE.*/
1810 /**/
1811 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
1812 /*    it under the terms of the GNU General Public License as published by*/
1813 /*    the Free Software Foundation; either version 2 of the License, or*/
1814 /*    (at your option) any later version.*/
1815 /**/
1816 /*    T-COFFEE is distributed in the hope that it will be useful,*/
1817 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
1818 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
1819 /*    GNU General Public License for more details.*/
1820 /**/
1821 /*    You should have received a copy of the GNU General Public License*/
1822 /*    along with Foobar; if not, write to the Free Software*/
1823 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
1824 /*...............................................                                                                                      |*/
1825 /*  If you need some more information*/
1826 /*  cedric.notredame@europe.com*/
1827 /*...............................................                                                                                                                                     |*/
1828 /**/
1829 /**/
1830 /*      */
1831 /******************************COPYRIGHT NOTICE*******************************/