binaries/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( short *table, int *pointt );
  15 void makecompositiontable( short *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( short *table, int *pointt )
  25 {
  26         int value = 0;
  27         short tmp;
  28         int point;
  29         static short *memo = NULL;
  30         static int *ct = NULL;
  31         static int *cp;
  32
  33         if( !memo )
  34         {
  35                 memo = vcalloc( tsize+1, sizeof( short ) );
  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 tupble exist
  61 */
  62 void makecompositiontable( short *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       short *table1;
 195       table1=vcalloc ( tsize,sizeof (short));
 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->L)
 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, CL->el_size);
 479     CL=index_constraint_list (CL);
 480
 481     l1=strlen (A->seq_al[l_s[0][0]]);
 482     l2=strlen (A->seq_al[l_s[1][0]]);
 483
 484     n_diag=l1+l2-1;
 485     diag=declare_int ( n_diag+2, 3);
 486     for ( a=0; a<= n_diag; a++)diag[a][0]=a;
 487
 488     A->S=CL->S;
 489     code=seq2aln_pos (A, ns, l_s);
 490     pos =aln2pos_simple ( A,-1, ns, l_s);
 491
 492
 493     for (a=0; a<ns[0]; a++)
 494
 495         {
 496         s1=A->order[l_s[0][a]][0];
 497         for (b=0; b<ns[1]; b++)
 498             {
 499             s2=A->order[l_s[1][b]][0];
 500             for (c=CL->start_index[s1][s2], d=0; c<CL->end_index[s1][s2];c++, d++)
 501                 {
 502                 entry=extract_entry ( entry,c, CL);
 503                 if (s1==entry[SEQ1])
 504                     {
 505                     r1=code [s1][entry[R1]];
 506                     r2=code [s2][entry[R2]];
 507                     }
 508                 else if ( s2==entry[SEQ1])
 509                     {
 510                     r2=code [s2][entry[R1]];
 511                     r1=code [s1][entry[R2]];
 512                     }
 513
 514
 515                 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);
 516
 517                 }
 518             }
 519         }
 520
 521
 522     sort_int (diag+1, 2, 1,0, n_diag-1);
 523
 524     free_int (code,-1);
 525     free_int (pos, -1);
 526     return diag;
 527     }
 528
 529 int * flag_diagonals (int l1, int l2, int **sorted_diag, float T, int window)
 530     {
 531     int a, b, up, low,current_diag,n_diag;
 532     int * slopes;
 533     int *diag_list;
 534     double mean;
 535     double sd;
 536     int use_z_score=1;
 537
 538
 539     n_diag=l1+l2-1;
 540     mean=return_mean_int ( sorted_diag, n_diag+1, 1);
 541
 542     sd  =return_sd_int ( sorted_diag, n_diag+1, 1, (int)mean);
 543
 544     if ( T==0)
 545       {
 546       use_z_score=1;
 547       T=(((double)sorted_diag[n_diag][1]-mean)/sd)/25;
 548       }
 549
 550
 551     diag_list=vcalloc (l1+l2+1, sizeof (int));
 552     slopes=vcalloc ( n_diag+1, sizeof (int));
 553
 554     for ( a=n_diag; a>0; a--)
 555             {
 556             current_diag=sorted_diag[a][0];
 557
 558
 559             if ( !use_z_score && sorted_diag[a][1]>T)
 560                {
 561                    up=MAX(1,current_diag-window);
 562                    low=MIN(n_diag, current_diag+window);
 563                    for ( b=up; b<=low; b++)slopes[b]=1;
 564                }
 565             else if (use_z_score && ((double)sorted_diag[a][1]-mean)/sd>T)
 566               {
 567                 up=MAX(1,current_diag-window);
 568                 low=MIN(n_diag, current_diag+window);
 569                 for ( b=up; b<=low; b++)slopes[b]=1;
 570               }
 571             else break;
 572             }
 573
 574     for ( a=1, b=0; a<=n_diag; a++)
 575         {
 576             b+=slopes[a];
 577         }
 578
 579     slopes[1]=1;
 580     slopes[l1+l2-1]=1;
 581     slopes[l2]=1;
 582     for (a=0; a<= (l1+l2-1); a++)
 583         if ( slopes[a]){diag_list[++diag_list[0]]=a;}
 584
 585     vfree (slopes);
 586
 587     return diag_list;
 588     }
 589 int * extract_N_diag (int l1, int l2, int **sorted_diag, int n_chosen_diag, int window)
 590     {
 591     int a, b, up, low,current_diag,n_diag;
 592     int * slopes;
 593     int *diag_list;
 594
 595
 596     n_diag=l1+l2-1;
 597
 598     diag_list=vcalloc (l1+l2+1, sizeof (int));
 599     slopes=vcalloc ( n_diag+1, sizeof (int));
 600
 601
 602     for ( a=n_diag; a>0 && a>(n_diag-n_chosen_diag); a--)
 603             {
 604             current_diag=sorted_diag[a][0];
 605             up=MAX(1,current_diag-window);
 606             low=MIN(n_diag, current_diag+window);
 607
 608             for ( b=up; b<=low; b++)slopes[b]=1;
 609             }
 610
 611     /*flag bottom right*/
 612     up=MAX(1,1-window);low=MIN(n_diag,1+window);
 613     for ( a=up; a<=low; a++) slopes[a]=1;
 614
 615     /*flag top left */
 616     up=MAX(1,(l1+l2-1)-window);low=MIN(n_diag,(l1+l2-1)+window);
 617     for ( a=up; a<=low; a++) slopes[a]=1;
 618
 619
 620     /*flag MAIN DIAG SEQ1*/
 621     up=MAX(1,l1-window);low=MIN(n_diag,l1+window);
 622     for ( a=up; a<=low; a++) slopes[a]=1;
 623
 624     /*flag MAIN DIAG SEQ2*/
 625     up=MAX(1,l2-window);low=MIN(n_diag,l2+window);
 626     for ( a=up; a<=low; a++) slopes[a]=1;
 627
 628
 629     for (a=0; a<= (l1+l2-1); a++)
 630         if ( slopes[a]){diag_list[++diag_list[0]]=a;}
 631
 632     vfree (slopes);
 633     return diag_list;
 634     }
 635
 636
 637
 638
 639 int cfasta_gotoh_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
 640     {
 641 /*TREATMENT OF THE TERMINAL GAP PENALTIES*/
 642 /*TG_MODE=0---> gop and gep*/
 643 /*TG_MODE=1---> ---     gep*/
 644 /*TG_MODE=2---> ---     ---*/
 645
 646
 647         int maximise;
 648
 649 /*VARIABLES FOR THE MULTIPLE SEQUENCE ALIGNMENT*/
 650         int **tot_diag;
 651
 652         int *diag;
 653         int ktup;
 654         static int n_groups;
 655         static char **group_list;
 656         int score, new_score;
 657         int n_chosen_diag=20;
 658         int step;
 659         int max_n_chosen_diag;
 660         int l1, l2;
 661         /********Prepare Penalties******/
 662
 663
 664         maximise=CL->maximise;
 665         ktup=CL->ktup;
 666
 667         /********************************/
 668
 669
 670
 671
 672         if ( !group_list)
 673            {
 674
 675                group_list=make_group_aa (&n_groups, CL->matrix_for_aa_group);
 676            }
 677
 678         l1=strlen (A->seq_al[l_s[0][0]]);
 679         l2=strlen (A->seq_al[l_s[1][0]]);
 680
 681         if ( !CL->fasta_step)
 682             {
 683             step=MIN(l1,l2);
 684             step=(int) log ((double)MAX(step, 1));
 685             step=MAX(step, 20);
 686             }
 687         else
 688             {
 689                 step=CL->fasta_step;
 690             }
 691
 692
 693         tot_diag=evaluate_diagonals ( A, ns, l_s, CL, maximise,n_groups,group_list, ktup);
 694
 695
 696         max_n_chosen_diag=strlen (A->seq_al[l_s[0][0]])+strlen (A->seq_al[l_s[1][0]])-1;
 697         /*max_n_chosen_diag=(int)log10((double)(l1+l2))*10;*/
 698
 699         n_chosen_diag+=step;
 700         n_chosen_diag=MIN(n_chosen_diag, max_n_chosen_diag);
 701
 702
 703         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);
 704
 705
 706         score    =make_fasta_gotoh_pair_wise ( A, ns, l_s, CL, diag);
 707
 708         new_score=0;
 709         vfree ( diag);
 710
 711
 712         while (new_score!=score && n_chosen_diag< max_n_chosen_diag    )
 713           {
 714
 715
 716             score=new_score;
 717
 718             ungap_sub_aln ( A, ns[0], l_s[0]);
 719             ungap_sub_aln ( A, ns[1], l_s[1]);
 720
 721
 722             n_chosen_diag+=step;
 723             n_chosen_diag=MIN(n_chosen_diag, max_n_chosen_diag);
 724
 725
 726             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);
 727             new_score=make_fasta_gotoh_pair_wise (  A, ns, l_s, CL, diag);
 728
 729             vfree ( diag);
 730
 731           }
 732
 733         score=new_score;
 734         free_int (tot_diag, -1);
 735
 736         return score;
 737     }
 738
 739 int fasta_gotoh_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
 740     {
 741 /*TREATMENT OF THE TERMINAL GAP PENALTIES*/
 742 /*TG_MODE=0---> gop and gep*/
 743 /*TG_MODE=1---> ---     gep*/
 744 /*TG_MODE=2---> ---     ---*/
 745
 746
 747         int maximise;
 748
 749 /*VARIABLES FOR THE MULTIPLE SEQUENCE ALIGNMENT*/
 750         int **tot_diag;
 751         int *diag;
 752         int ktup;
 753         float diagonal_threshold;
 754         static int n_groups;
 755         static char **group_list;
 756         int score;
 757         /********Prepare Penalties******/
 758
 759
 760         maximise=CL->maximise;
 761         ktup=CL->ktup;
 762         diagonal_threshold=CL->diagonal_threshold;
 763         /********************************/
 764
 765
 766
 767         if ( !group_list)
 768            {
 769                group_list=make_group_aa (&n_groups, CL->matrix_for_aa_group);
 770            }
 771
 772
 773         tot_diag=evaluate_diagonals ( A, ns, l_s, CL, maximise,n_groups,group_list, ktup);
 774
 775         if (  !CL->fasta_step)
 776           {
 777             diag=flag_diagonals (strlen(A->seq_al[l_s[0][0]]),strlen(A->seq_al[l_s[1][0]]), tot_diag,diagonal_threshold,0);
 778           }
 779
 780         else
 781           {
 782
 783             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);
 784
 785           }
 786         score=make_fasta_gotoh_pair_wise ( A, ns, l_s, CL, diag);
 787
 788         free_int (tot_diag, -1);
 789         vfree (diag);
 790         return score;
 791     }
 792 int very_fast_gotoh_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
 793     {
 794 /*TREATMENT OF THE TERMINAL GAP PENALTIES*/
 795 /*TG_MODE=0---> gop and gep*/
 796 /*TG_MODE=1---> ---     gep*/
 797 /*TG_MODE=2---> ---     ---*/
 798
 799
 800         int maximise;
 801 /*VARIABLES FOR THE MULTIPLE SEQUENCE ALIGNMENT*/
 802         int **tot_diag;
 803         int *diag;
 804         int ktup;
 805         static int n_groups;
 806         static char **group_list;
 807         int score;
 808         /********Prepare Penalties******/
 809
 810
 811         maximise=CL->maximise;
 812         ktup=CL->ktup;
 813         /********************************/
 814
 815
 816         if ( !group_list)
 817            {
 818
 819                group_list=make_group_aa (&n_groups, CL->matrix_for_aa_group);
 820            }
 821
 822         CL->use_fragments=0;
 823         tot_diag=evaluate_diagonals ( A, ns, l_s, CL, maximise,n_groups,group_list, ktup);
 824
 825         /*Note: 20 diagonals. 5 shadows on each side: tunned on Hom39, 2/2/04 */
 826         diag=extract_N_diag (strlen (A->seq_al[l_s[0][0]]),strlen (A->seq_al[l_s[1][0]]), tot_diag,20,5);
 827         score=make_fasta_gotoh_pair_wise ( A, ns, l_s, CL, diag);
 828         free_int (tot_diag, -1);
 829         vfree (diag);
 830         return score;
 831     }
 832 int make_fasta_gotoh_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL, int *diag)
 833     {
 834 /*TREATMENT OF THE TERMINAL GAP PENALTIES*/
 835 /*TG_MODE=0---> gop and gep*/
 836 /*TG_MODE=1---> ---     gep*/
 837       /*TG_MODE=2---> ---     ---*/
 838
 839
 840         int TG_MODE, gop, l_gop, gep,l_gep, maximise;
 841
 842 /*VARIABLES FOR THE MULTIPLE SEQUENCE ALIGNMENT*/
 843         int a, b,c,k, t;
 844         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;
 845         int su, in, de, tr;
 846
 847         int **C, **D, **I, **trace, **pos0, **LD;
 848         int lenal[2], len;
 849         char *buffer, *char_buf;
 850         char **aln, **al;
 851
 852         /********Prepare Penalties******/
 853         gop=CL->gop*SCORE_K;
 854         gep=CL->gep*SCORE_K;
 855         TG_MODE=CL->TG_MODE;
 856         maximise=CL->maximise;
 857
 858
 859         /********************************/
 860
 861
 862         n_diag=diag[0];
 863
 864
 865
 866        l1=lenal[0]=strlen (A->seq_al[l_s[0][0]]);
 867        l2=lenal[1]=strlen (A->seq_al[l_s[1][0]]);
 868
 869        if ( getenv ("DEBUG_TCOFFEE"))fprintf ( stderr, "\n\tNdiag=%d%%  ", (diag[0]*100)/(l1+l2));
 870
 871         /*diag:
 872           diag[1..n_diag]--> flaged diagonal in order;
 873           diag[0]=0--> first diagonal;
 874           diag[n_diag+1]=l1+l2-1;
 875         */
 876
 877         /*numeration of the diagonals strats from the bottom right [1...l1+l2-1]*/
 878         /*sequence s1 is vertical and seq s2 is horizontal*/
 879         /*D contains the best Deletion  in S2==>comes from diagonal N+1*/
 880         /*I contains the best insertion in S2=> comes from diagonal N-1*/
 881
 882
 883
 884
 885
 886        C=declare_int (lenal[0]+lenal[1]+1, n_diag+2);
 887        D=declare_int (lenal[0]+lenal[1]+1, n_diag+2);
 888        LD=declare_int (lenal[0]+lenal[1]+1, n_diag+2);
 889        I=declare_int (lenal[0]+lenal[1]+1, n_diag+2);
 890        trace=declare_int (lenal[0]+lenal[1]+1, n_diag+2);
 891
 892
 893        al=declare_char (2,lenal[0]+lenal[1]+lenal[1]+1);
 894
 895        len= MAX(lenal[0],lenal[1])+1;
 896        buffer=vcalloc ( 2*len, sizeof (char));
 897        char_buf= vcalloc (2*len, sizeof (char));
 898
 899        pos0=aln2pos_simple ( A,-1, ns, l_s);
 900        C[0][0]=0;
 901
 902        t=(TG_MODE==0)?gop:0;
 903        for ( j=1; j<= n_diag; j++)
 904             {
 905                 l_gop=(TG_MODE==0)?gop:0;
 906                 l_gep=(TG_MODE==2)?0:gep;
 907
 908
 909
 910                 if ( (diag[j]-lenal[0])<0 )
 911                     {
 912                     trace[0][j]=UNDEFINED;
 913                     continue;
 914                     }
 915                 C[0][j]=(diag[j]-lenal[0])*l_gep +l_gop;
 916                 D[0][j]=(diag[j]-lenal[0])*l_gep +l_gop+gop;
 917             }
 918        D[0][j]=D[0][j-1]+gep;
 919
 920
 921        t=(TG_MODE==0)?gop:0;
 922        for ( i=1; i<=lenal[0]; i++)
 923            {
 924                 l_gop=(TG_MODE==0)?gop:0;
 925                 l_gep=(TG_MODE==2)?0:gep;
 926
 927                 C[i][0]=C[i][n_diag+1]=t=t+l_gep;
 928                 I[i][0]=D[i][n_diag+1]=t+    gop;
 929
 930                 for ( j=1; j<=n_diag; j++)
 931                     {
 932                         C[i][j]=C[i][0];
 933                         D[i][j]=I[i][j]=I[i][0];
 934                     }
 935
 936                 for (eg=0, j=1; j<=n_diag; j++)
 937                     {
 938
 939                         pos_j=diag[j]-lenal[0]+i;
 940                         if (pos_j<=0 || pos_j>l2 )
 941                             {
 942                             trace[i][j]=UNDEFINED;
 943                             continue;
 944                             }
 945                         sub=(CL->get_dp_cost) ( A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],pos_j-1, CL );
 946
 947                     /*1 identify the best insertion in S2:*/
 948                         l_gop=(i==lenal[0])?((TG_MODE==0)?gop:0):gop;
 949                         l_gep=(i==lenal[0])?((TG_MODE==2)?0:gep):gep;
 950                         len=(j==1)?0:(diag[j]-diag[j-1]);
 951                         if ( a_better_than_b(I[i][j-1], C[i][j-1]+l_gop, maximise))eg++;
 952                         else eg=1;
 953                         I[i][j]=best_of_a_b (I[i][j-1], C[i][j-1]+l_gop, maximise)+len*l_gep;
 954
 955                     /*2 Identify the best deletion in S2*/
 956                         l_gop=(pos_j==lenal[1])?((TG_MODE==0)?gop:0):gop;
 957                         l_gep=(pos_j==lenal[1])?((TG_MODE==2)?0:gep):gep;
 958
 959                         len=(j==n_diag)?0:(diag[j+1]-diag[j]);
 960                         delta_i=((i-len)>0)?(i-len):0;
 961
 962                         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;}
 963                         else {LD[i][j]=1;}
 964                         D[i][j]=best_of_a_b (D[delta_i][j+1],C[delta_i][j+1]+l_gop, maximise)+len*l_gep;
 965
 966
 967                         /*Identify the best way*/
 968                         /*
 969                         score=C[i][j]=best_int ( 3, maximise, &fop, I[i][j], C[i-1][j]+sub, D[i][j]);
 970                         fop-=1;
 971                         if ( fop<0)trace[i][j]=fop*eg;
 972                         else if ( fop>0 ) {trace[i][j]=fop*LD[i][j];}
 973                         else if ( fop==0) trace[i][j]=0;
 974                         */
 975
 976                         su=C[i-1][j]+sub;
 977                         in=I[i][j];
 978                         de=D[i][j];
 979
 980                         /*HERE ("%d %d %d", su, in, de);*/
 981                         if (su>=in && su>=de)
 982                           {
 983                             score=su;
 984                             tr=0;
 985                           }
 986                         else if (in>=de)
 987                           {
 988                             score=in;
 989                             tr=-eg;
 990                           }
 991                         else
 992                           {
 993                             score=de;
 994                             tr=LD[i][j];
 995                           }
 996                         trace[i][j]=tr;
 997                         C[i][j]=score;
 998
 999
1000                         last_i=i;
1001                         last_j=j;
1002                     }
1003             }
1004
1005
1006        /*
1007                     [0][Positive]
1008                      ^     ^
1009                      |    /
1010                      |   /
1011                      |  /
1012                      | /
1013                      |/
1014        [Neg]<-------[*]
1015         */
1016
1017
1018         i=last_i;
1019         j=last_j;
1020
1021
1022
1023         ala=alb=0;
1024         match1=match2=0;
1025         while (!(match1==l1 && match2==l2))
1026               {
1027
1028
1029                   if ( match1==l1)
1030                      {
1031                          len=l2-match2;
1032                          for ( a=0; a< len; a++)
1033                              {
1034                              al[0][ala++]=0;
1035                              al[1][alb++]=1;
1036                              match2++;
1037                              }
1038                          k=0;
1039                          break;
1040
1041                          /*k=-(j-1);*/
1042
1043                      }
1044                   else if ( match2==l2)
1045                      {
1046                          len=l1-match1;
1047                          for ( a=0; a< len; a++)
1048                              {
1049                              al[0][ala++]=1;
1050                              al[1][alb++]=0;
1051                              match1++;
1052                              }
1053                          k=0;
1054                          break;
1055                          /*k= n_diag-j;*/
1056                      }
1057                   else
1058                       {
1059                           k=trace[i][j];
1060                       }
1061
1062
1063                   if ( k==0)
1064                              {
1065                                  if ( match2==l2 || match1==l1);
1066                                  else
1067                                     {
1068
1069                                     al[0][ala++]=1;
1070                                     al[1][alb++]=1;
1071                                     i--;
1072                                     match1++;
1073                                     match2++;
1074                                     }
1075                              }
1076                   else if ( k>0)
1077                              {
1078
1079                              len=diag[j+k]-diag[j];
1080                              for ( a=0; a<len; a++)
1081                                  {
1082                                      if ( match1==l1)break;
1083                                      al[0][ala++]=1;
1084                                      al[1][alb++]=0;
1085                                      match1++;
1086                                  }
1087                              i-=len;
1088                              j+=k;
1089                              }
1090                   else if ( k<0)
1091                              {
1092                              k*=-1;
1093                              len=diag[j]-diag[j-k];
1094                              for ( a=0; a<len; a++)
1095                                  {
1096                                      if ( match2==l2)break;
1097                                      al[0][ala++]=0;
1098                                      al[1][alb++]=1;
1099                                      match2++;
1100                                  }
1101
1102
1103                              j-=k;
1104                              }
1105               }
1106
1107         LEN=ala;
1108         c=LEN-1;
1109         invert_list_char ( al[0], LEN);
1110         invert_list_char ( al[1], LEN);
1111         if ( A->declared_len<=LEN)A=realloc_aln2  ( A,A->max_n_seq, 2*LEN);
1112         aln=A->seq_al;
1113
1114         for ( c=0; c< 2; c++)
1115             {
1116             for ( a=0; a< ns[c]; a++)
1117                 {
1118                 ch=0;
1119                 for ( b=0; b< LEN; b++)
1120                     {
1121                     if (al[c][b]==1)
1122                         char_buf[b]=aln[l_s[c][a]][ch++];
1123                     else
1124                         char_buf[b]='-';
1125                    }
1126                 char_buf[b]='\0';
1127                 sprintf (aln[l_s[c][a]],"%s", char_buf);
1128                 }
1129              }
1130
1131
1132         A->len_aln=LEN;
1133         A->nseq=ns[0]+ns[1];
1134
1135         free_int (pos0, -1);
1136         free_int (C, -1);
1137         free_int (D, -1);
1138         free_int (I, -1);
1139         free_int (trace, -1);
1140         free_int (LD, -1);
1141         free_char ( al, -1);
1142         vfree(buffer);
1143         vfree(char_buf);
1144
1145
1146         return score;
1147     }
1148
1149 int hasch_seq(char *seq, int **hs, int **lu,int ktup,char *alp)
1150     {
1151         static int a[10];
1152
1153         int i,j,l,limit,code,flag;
1154         char residue;
1155
1156         int alp_lu[10000];
1157         int alp_size;
1158
1159         alp_size=alp[0];
1160         alp++;
1161
1162
1163
1164         for ( i=0; i< alp_size; i++)
1165             {
1166               alp_lu[(int)alp[i]]=i;
1167             }
1168
1169
1170
1171         l=strlen (seq);
1172         limit = (int)   pow((double)(alp_size+1),(double)ktup);
1173         hs[0]=vcalloc ( l+1,sizeof (int));
1174         lu[0]=vcalloc ( limit+1, sizeof(int));
1175
1176
1177         if ( l==0)myexit(EXIT_FAILURE);
1178
1179         for (i=1;i<=ktup;i++)
1180            a[i] = (int) pow((double)(alp_size+1),(double)(i-1));
1181
1182
1183         for(i=1;i<=(l-ktup+1);++i)
1184                 {
1185                 code=0;
1186                 flag=FALSE;
1187                 for(j=1;j<=ktup;++j)
1188                    {
1189                    if (is_gap(seq[i+j-2])){flag=TRUE;break;}
1190                    else residue=alp_lu[(int)seq[i+j-2]];
1191                    code+=residue*a[j];
1192                    }
1193
1194                 if ( flag)continue;
1195                 ++code;
1196
1197                 if (lu[0][code])hs[0][i]=lu[0][code];
1198                 lu[0][code]=i;
1199                 }
1200         return 0;
1201     }
1202
1203
1204
1205 /*********************************************************************/
1206 /*                                                                   */
1207 /*                         KTUP_DP                                   */
1208 /*                                                                   */
1209 /*                                                                   */
1210 /*********************************************************************/
1211
1212 /**************Hasch DAta Handling*******************************************************/
1213
1214 struct Hasch_data * free_ktup_hasch_data (struct Hasch_data *d);
1215 struct Hasch_data * declare_ktup_hasch_data (struct Hasch_entry *e);
1216 struct Hasch_data * allocate_ktup_hasch_data (struct Hasch_data *e, int action);
1217
1218 struct Hasch_data
1219 {
1220  int *list;
1221 };
1222 typedef struct Hasch_data Hasch_data;
1223 struct Hasch_data * free_ktup_hasch_data (struct Hasch_data *d)
1224 {
1225   return allocate_ktup_hasch_data (d, FREE);
1226 }
1227 struct Hasch_data * declare_ktup_hasch_data (struct Hasch_entry *e)
1228 {
1229   e->data=allocate_ktup_hasch_data (NULL,DECLARE);
1230   return e->data;
1231 }
1232
1233 struct Hasch_data * allocate_ktup_hasch_data (struct Hasch_data *e, int action)
1234 {
1235   static struct Hasch_data **heap;
1236   static int heap_size, free_heap, a;
1237
1238   if ( action == 100)
1239     {
1240       fprintf ( stderr, "\nHeap size: %d, Free Heap: %d", heap_size, free_heap);
1241       return NULL;
1242     }
1243   else if ( action==DECLARE)
1244     {
1245       if ( free_heap==0)
1246         {
1247           free_heap=100;
1248           heap_size+=free_heap;
1249           heap=vrealloc (heap,heap_size*sizeof (struct Hasch_entry *));
1250           for ( a=0; a<free_heap; a++)
1251             {
1252               (heap[a])=vcalloc ( 1, sizeof ( struct Hasch_entry *));
1253               (heap[a])->list=vcalloc ( 10, sizeof (int));
1254               (heap[a])->list[0]=10;
1255             }
1256         }
1257       return heap[--free_heap];
1258     }
1259   else if ( action==FREE)
1260     {
1261       heap[free_heap++]=e;
1262       e->list[1]=0;
1263       return NULL;
1264     }
1265   return NULL;
1266 }
1267
1268
1269 /**************Hasch DAta Handling*******************************************************/
1270
1271 int precomputed_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
1272     {
1273       int l1, l2, a, b, c;
1274       int nid=0, npos=0, id;
1275       int r1, r2, s1, s2;
1276
1277       l1=strlen(A->seq_al[l_s[0][0]]);
1278       l2=strlen(A->seq_al[l_s[1][0]]);
1279       if (l1!=l2)
1280         {
1281           fprintf ( stderr, "\nERROR: improper use of the function precomputed pairwise:[FATAL:%s]", PROGRAM);
1282           crash ("");
1283         }
1284       else if ( l1==0)
1285         {
1286           A->score_aln=A->score=0;
1287           return 0;
1288         }
1289
1290       for (npos=0, nid=0, a=0; a< ns[0]; a++)
1291         {
1292           s1=l_s[0][a];
1293
1294           for (b=0; b< ns[1]; b++)
1295             {
1296               s2=l_s[1][b];
1297               for ( c=0; c<l1; c++)
1298                 {
1299                 r1=A->seq_al[s1][c];
1300                 r2=A->seq_al[s2][c];
1301                 if ( is_gap(r1) || is_gap(r2));
1302                 else
1303                   {
1304                     npos++;
1305                     nid+=(r1==r2);
1306                   }
1307                 }
1308             }
1309         }
1310       id=(npos==0)?0:((nid*100)/npos);
1311       A->score=A->score_aln=id;
1312       return A->score;
1313     }
1314 int ktup_comparison_str ( char *seq1, char *seq2, const int ktup);
1315 int ktup_comparison_hasch ( char *i_seq1, char *i_seq2, const int ktup);
1316 int ktup_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
1317     {
1318       static char **gl;
1319       static int ng;
1320       char *seq1;
1321       char *seq2;
1322
1323       int min_len=10;
1324
1325
1326
1327       if ( !gl)
1328         gl=make_group_aa (&ng, "vasiliky");
1329
1330
1331       if ( ns[0]>1)seq1=sub_aln2cons_seq_mat (A, ns[0], l_s[0],"blosum62mt");
1332       else
1333         {
1334           seq1=vcalloc ( strlen (A->seq_al[l_s[0][0]])+1, sizeof (char));
1335           sprintf ( seq1, "%s",A->seq_al[l_s[0][0]]);
1336         }
1337       if ( ns[1]>1)seq2=sub_aln2cons_seq_mat (A, ns[1], l_s[1],"blosum62mt");
1338       else
1339         {
1340           seq2=vcalloc ( strlen (A->seq_al[l_s[1][0]])+1, sizeof (char));
1341           sprintf ( seq2, "%s",A->seq_al[l_s[1][0]]);
1342         }
1343
1344       if ( strlen (seq1)<min_len || strlen (seq2)<min_len)
1345         {
1346           Alignment *B;
1347
1348           ungap(seq1); ungap(seq2);
1349           B=align_two_sequences ( seq1, seq2, "blosum62mt",-10, -1, "myers_miller_pair_wise");
1350           A->score=A->score_aln=aln2sim(B, "idmat");
1351           free_aln (B);
1352           return A->score;
1353         }
1354       else
1355         {
1356
1357           string_convert (seq1, ng, gl);
1358           string_convert (seq2, ng, gl);
1359           A->score=A->score_aln=ktup_comparison (seq1,seq2, CL->ktup);
1360         }
1361
1362       vfree (seq1); vfree (seq2);
1363       return A->score;
1364     }
1365 int ktup_comparison( char *seq2, char *seq1, const int ktup)
1366 {
1367   return ktup_comparison_hasch ( seq2, seq1, ktup);
1368 }
1369 int ktup_comparison_str ( char *seq2, char *seq1, const int ktup)
1370 {
1371   int a,l1, l2,c1, c2, end, start;
1372   char *s1, *s2;
1373   double score=0;
1374   int max_dist=-1;
1375
1376   if ( max_dist==-1)max_dist=MAX((strlen (seq1)),(strlen (seq2)));
1377   l1=strlen (seq1)-ktup;
1378   l2=strlen (seq2);
1379
1380
1381   for ( a=0; a< l1; a++)
1382     {
1383       c1=seq1[a+ktup];seq1[a+ktup]='\0';
1384       s1=seq1+a;
1385
1386       start=((a-max_dist)<0)?0:a-max_dist;
1387       end=((a+max_dist)>=l2)?l2:a+max_dist;
1388
1389       c2=seq2[end];seq2[end]='\0';
1390       s2=seq2+start;
1391
1392       score+=(strstr(s2, s1)!=NULL)?1:0;
1393
1394       seq1[a+ktup]=c1;
1395       seq2[end]=c2;
1396     }
1397   score/=(l1==0)?1:l1;
1398   score=((log(0.1+score)-log(0.1))/(log(1.1)-log(0.1)));
1399
1400   return score*100;
1401
1402 }
1403 int ktup_comparison_hasch ( char *i_seq1, char *i_seq2, const int ktup)
1404 {
1405   /*Ktup comparison adapted from Rob Edgar, NAR, vol32, No1, 381, 2004*/
1406   /*1: hasch sequence 1
1407     2: Count the number of seq2 ktup found in seq1
1408   */
1409
1410   char c;
1411   int key;
1412
1413   static HaschT*H1;
1414   static char *pseq;
1415   Hasch_entry *e;
1416   char *s;
1417   int l, ls;
1418   int p, a, max_dist=-1;
1419   double score=0;
1420
1421
1422
1423   if (!strm (i_seq1, pseq))
1424     {
1425       if (H1)
1426         {
1427           hdestroy (H1, declare_ktup_hasch_data, free_ktup_hasch_data);
1428           string2key (NULL, NULL);
1429         }
1430       H1=hasch_sequence ( i_seq1, ktup);
1431       vfree (pseq);pseq=vcalloc ( strlen (i_seq1)+1, sizeof (char));
1432       sprintf ( pseq, "%s", i_seq1);
1433     }
1434
1435   ls=l=strlen (i_seq2);
1436   s=i_seq2;
1437   p=0;
1438   while (ls>ktup)
1439     {
1440       c=s[ktup];s[ktup]='\0';
1441       key=string2key (s, NULL);
1442       e=hsearch (H1,key,FIND, declare_ktup_hasch_data, free_ktup_hasch_data);
1443
1444       if ( e==NULL);
1445       else if ( max_dist==-1)score++;
1446       else
1447         {
1448           for ( a=1; a<=(e->data)->list[1]; a++)
1449             if (FABS((p-(e->data)->list[a]))<=max_dist)
1450               {score++; break;}
1451         }
1452       s[ktup]=c;s++;p++;ls--;
1453     }
1454   score/=(l-ktup);
1455   score=(log(0.1+score)-log(0.1))/(log(1.1)-log(0.1));
1456
1457   if ( score>100) score=100;
1458   return (int)(score*100);
1459 }
1460
1461 HaschT* hasch_sequence ( char *seq1, int ktup)
1462 {
1463   char c;
1464   int key, offset=0, ls;
1465   HaschT *H;
1466   Hasch_entry *e;
1467
1468   H=hcreate ( strlen (seq1), declare_ktup_hasch_data, free_ktup_hasch_data);
1469   ls=strlen (seq1);
1470   while (ls>=(ktup))
1471     {
1472       c=seq1[ktup];seq1[ktup]='\0';
1473       key=string2key (seq1, NULL);
1474       e=hsearch (H,key,FIND, declare_ktup_hasch_data, free_ktup_hasch_data);
1475
1476       if (e==NULL)
1477         {
1478          e=hsearch (H,key,ADD,declare_ktup_hasch_data,free_ktup_hasch_data);
1479          (e->data)->list[++(e->data)->list[1]+1]=offset;
1480         }
1481       else
1482         {
1483           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));}
1484           (e->data)->list[++(e->data)->list[1]+1]=offset;
1485         }
1486        seq1[ktup]=c;seq1++;ls--;
1487        offset++;
1488     }
1489   return H;
1490 }
1491
1492
1493
1494 char *dayhoff_translate (char *seq1)
1495 {
1496 int l, a, c;
1497 l=strlen (seq1);
1498  for ( a=0; a< l; a++)
1499   {
1500     c=tolower(seq1[a]);
1501     if ( strchr ("agpst", c))seq1[a]='a';
1502     else if (strchr ("denq", c))seq1[a]='d';
1503     else if (strchr ("fwy", c))seq1[a]='f';
1504     else if (strchr ("hkr", c))seq1[a]='h';
1505     else if (strchr ("ilmv", c))seq1[a]='i';
1506   }
1507 return seq1;
1508 }
1509
1510 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)
1511 {
1512   /*Ktup comparison as in Rob Edgar, NAR, vol32, No1, 381, 2004*/
1513   char character;
1514   int key,ls;
1515   HaschT*H1, *H2;
1516   Hasch_entry *e1, *e2;
1517   char *s, *sb, *seq1, *seq2;
1518   int l1, l2;
1519   int score=0;
1520   int **diag,n_diag, ktup1, ktup2,a,b,c,d, **pos;
1521   int n_dots=0;
1522
1523   pos=aln2pos_simple ( A,-1, ns, l_s);
1524
1525   seq1=aln2cons_maj (A, ns[0], l_s[0], n_groups, group_list);
1526   seq2=aln2cons_maj (A, ns[1], l_s[1], n_groups, group_list);
1527   l1=strlen (seq1);
1528   l2=strlen (seq2);
1529   n_diag=l1+l2-1;
1530
1531
1532   diag=declare_int (n_diag+2, 3);
1533   for ( a=0; a<n_diag+2; a++)diag[a][0]=a;
1534
1535   H1=hasch_sequence ( seq1, ktup);
1536   H2=hasch_sequence ( seq2, ktup);
1537   s=sb=vcalloc (strlen (seq1)+strlen (seq2)+1, sizeof (char));
1538   sprintf (s, "%s%s", seq1, seq2);
1539
1540   ls=strlen(s);
1541   while (ls>=(ktup))
1542     {
1543       character=s[ktup];s[ktup]='\0';
1544       key=string2key (s, NULL);
1545       e1=hsearch (H1,key,FIND,declare_ktup_hasch_data, free_ktup_hasch_data);
1546       e2=hsearch (H2,key,FIND,declare_ktup_hasch_data, free_ktup_hasch_data);
1547       if ( !e2 || !e1);
1548       else
1549         {
1550
1551           for (b=2; b<(e1->data)->list[1]+2; b++)
1552             for (c=2; c<(e2->data)->list[1]+2; c++)
1553               {
1554
1555                 ktup1=(e1->data)->list[b];
1556                 ktup2=(e2->data)->list[c];
1557                 diag[(ktup2-ktup1)+l1][2]++;
1558                 for (score=0, d=0; d<ktup; d++)
1559                   score+=(CL->get_dp_cost) ( A, pos, ns[0], l_s[0], ktup1+d, pos,ns[1], l_s[1], ktup2+d, CL);
1560                 diag[(ktup2-ktup1)+l1][1]+=score;
1561                 n_dots++;
1562               }
1563           (e1->data)->list[1]=(e2->data)->list[1]=0;
1564         }
1565       s[ktup]=character;s++;ls--;
1566     }
1567
1568   sort_int (diag+1, 2, 1,0,n_diag-1);
1569
1570   hdestroy (H1,declare_ktup_hasch_data, free_ktup_hasch_data); hdestroy (H2,declare_ktup_hasch_data, free_ktup_hasch_data);
1571   vfree (seq1); vfree (seq2);vfree (sb);free_int (pos, -1);
1572   return diag;
1573 }
1574  /*********************************************************************/
1575 /*                                                                   */
1576 /*                         OLD FUNCTIONS                              */
1577 /*                                                                   */
1578 /*                                                                   */
1579 /*********************************************************************/
1580 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)
1581     {
1582    /*
1583     Reads in an alignmnent A, with two groups of sequences marked.
1584     1-Turn each group into a conscensus, using the group list identifier.
1585                -if the group list is left empty original symbols are used
1586     2-hasch the two sequences
1587     3-score each diagonal, sort the list and return it (diag_list)
1588
1589         diag_list:
1590
1591    */
1592
1593     char *seq1, *seq2, *alphabet=NULL;
1594     int a,b,l1, l2, n_ktup,pos_ktup1, pos_ktup2, **pos;
1595     int *hasched_seq1, *hasched_seq2,*lu_seq1,*lu_seq2;
1596     int n_diag, **diag, current_diag, n_dots;
1597     static char *buf;
1598     pos=aln2pos_simple ( A,-1, ns, l_s);
1599
1600
1601     seq1=aln2cons_seq (A, ns[0], l_s[0], n_groups, group_list);
1602     seq2=aln2cons_seq (A, ns[1], l_s[1], n_groups, group_list);
1603
1604
1605
1606
1607     alphabet=get_alphabet (seq1,alphabet);
1608     alphabet=get_alphabet (seq2,alphabet);
1609
1610     l1=strlen ( seq1);
1611     l2=strlen ( seq2);
1612
1613     n_diag=l1+l2-1;
1614     diag=declare_int ( n_diag+2, 3);
1615     n_ktup=(int)pow ( (double)alphabet[0]+1, (double)ktup);
1616
1617
1618     hasch_seq(seq1, &hasched_seq1, &lu_seq1,ktup, alphabet);
1619     hasch_seq(seq2, &hasched_seq2, &lu_seq2,ktup, alphabet);
1620
1621
1622
1623
1624     /*EVALUATE THE DIAGONALS*/
1625     for ( a=0; a<= n_diag; a++)diag[a][0]=a;
1626     for ( n_dots=0,a=1; a<= n_ktup; a++)
1627         {
1628             pos_ktup1=lu_seq1[a];
1629             while (TRUE)
1630                   {
1631                   if (!pos_ktup1)break;
1632                   pos_ktup2=lu_seq2[a];
1633                   while (pos_ktup2)
1634                             {
1635                             current_diag=(pos_ktup2-pos_ktup1+l1);
1636                             for ( b=0; b< ktup; b++)
1637                                 {
1638                                     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);
1639                                     n_dots++;
1640
1641                                 }
1642                             diag[current_diag][2]++;
1643                             pos_ktup2=hasched_seq2[pos_ktup2];
1644                             }
1645                   pos_ktup1=hasched_seq1[pos_ktup1];
1646                   }
1647
1648         }
1649     if ( n_dots==0)
1650        {
1651            if ( !buf)
1652                {
1653                buf=vcalloc ( 30, sizeof (30));
1654                sprintf ( buf, "abcdefghijklmnopqrstuvwxyz");
1655                }
1656             vfree ( hasched_seq1);
1657             vfree ( hasched_seq2);
1658             vfree (lu_seq1);
1659             vfree (lu_seq2);
1660            return evaluate_diagonals_with_ktup ( A,ns,l_s, CL,maximise,1,&buf,1);
1661        }
1662
1663
1664     sort_int (diag+1, 2, 1,0, n_diag-1);
1665     vfree (seq1);
1666     vfree (seq2);
1667     vfree (alphabet);
1668     vfree ( hasched_seq1);
1669     vfree ( hasched_seq2);
1670     vfree (lu_seq1);
1671     vfree (lu_seq2);
1672     free_int (pos, -1);
1673     return diag;
1674     }
1675 /////////////////////////////////////////////////////////////////
1676
1677 Constraint_list * hasch2constraint_list (Sequence*S, Constraint_list *CL)
1678 {
1679   int a,b,c, n;
1680   SeqHasch h,*H=NULL;
1681   int *entry;
1682   int ktup=2;
1683
1684
1685   entry=vcalloc ( CL->entry_len, sizeof (int));
1686
1687   for (a=0; a<S->nseq; a++)
1688     {
1689       H=seq2hasch (a, S->seq[a],ktup,H);
1690     }
1691
1692   n=1;
1693   while (H[n])
1694     {
1695       h=H[n];
1696
1697       for (a=0; a<h->n-2; a+=2)
1698         {
1699           for (b=a+2; b<h->n; b+=2)
1700             {
1701
1702               if (h->l[a]==h->l[b])continue;
1703               else
1704                 {
1705                   for (c=0; c<ktup; c++)
1706                     {
1707                       entry[SEQ1]=h->l[a];
1708                       entry[SEQ2]=h->l[b];
1709                       entry[R1]=h->l[a+1]+c;
1710                       entry[R2]=h->l[b+1]+c;
1711                       entry[WE]=100;
1712                       add_entry2list (entry,CL);
1713                     }
1714                 }
1715             }
1716         }
1717       n++;
1718     }
1719
1720   return CL;
1721 }
1722 SeqHasch *cleanhasch       (SeqHasch *H)
1723 {
1724   int n=1;
1725   SeqHasch *N;
1726   N=vcalloc (2, sizeof (SeqHasch));
1727   N[0]=H[0];
1728
1729   while (H[n])
1730     {
1731       (H[n])->n=0;
1732       vfree ((H[n])->l);
1733       (H[n])->l=NULL;
1734       n++;
1735     }
1736   vfree (H);
1737   return N;
1738 }
1739 int hasch2sim        (SeqHasch *H, int nseq)
1740 {
1741   int n=1;
1742
1743   int a,cs, ps, ns;
1744   int id=0, tot=0;
1745
1746   while (H[n])
1747     {
1748       for (ps=-1,ns=0,a=0; a<(H[n])->n; a+=2)
1749         {
1750           //HERE ("%d--[%d %d]",n, (H[n])->l[a], (H[n])->l[a+1]);
1751           cs=(H[n])->l[a];
1752           if (cs!=ps)ns++;
1753           ps=cs;
1754         }
1755       n++;
1756       if (ns==nseq)id++;
1757       tot++;
1758     }
1759
1760   return (id*MAXID)/tot;
1761 }
1762 SeqHasch * seq2hasch (int i,char *seq, int ktup, SeqHasch *H)
1763 {
1764   int a,b,l, n=0;
1765   SeqHasch h;
1766
1767
1768   if (!H)
1769     {
1770       H=vcalloc (2, sizeof (SeqHasch));
1771       H[0]=vcalloc (1, sizeof (hseq));
1772       n=1;
1773     }
1774   else
1775     {
1776       n=0;
1777       while (H[++n]);
1778     }
1779
1780   l=strlen (seq);
1781   for (a=0; a<l-ktup; a++)
1782     {
1783       h=H[0];
1784       for (b=a; b<a+ktup; b++)
1785         {
1786           char r;
1787           r=seq[b];
1788           if (!h->hl[r])  h->hl[r]=vcalloc (1, sizeof (hseq));
1789           h=h->hl[r];
1790         }
1791       if (!h->l)
1792         {
1793
1794           h->n=2;
1795           h->l=vcalloc (2, sizeof (int));
1796           H=vrealloc (H,(n+2)*sizeof (SeqHasch));
1797           H[n]=h;
1798           n++;
1799         }
1800       else
1801         {
1802           h->n+=2;
1803           h->l=vrealloc (h->l, (h->n)*sizeof (int));
1804         }
1805
1806       h->l[h->n-2]=i;
1807       h->l[h->n-1]=a;
1808     }
1809   return H;
1810 }
1811
1812 /*********************************COPYRIGHT NOTICE**********************************/
1813 /*© Centro de Regulacio Genomica */
1814 /*and */
1815 /*Cedric Notredame */
1816 /*Tue Oct 27 10:12:26 WEST 2009. */
1817 /*All rights reserved.*/
1818 /*This file is part of T-COFFEE.*/
1819 /**/
1820 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
1821 /*    it under the terms of the GNU General Public License as published by*/
1822 /*    the Free Software Foundation; either version 2 of the License, or*/
1823 /*    (at your option) any later version.*/
1824 /**/
1825 /*    T-COFFEE is distributed in the hope that it will be useful,*/
1826 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
1827 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
1828 /*    GNU General Public License for more details.*/
1829 /**/
1830 /*    You should have received a copy of the GNU General Public License*/
1831 /*    along with Foobar; if not, write to the Free Software*/
1832 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
1833 /*...............................................                                                                                      |*/
1834 /*  If you need some more information*/
1835 /*  cedric.notredame@europe.com*/
1836 /*...............................................                                                                                                                                     |*/
1837 /**/
1838 /**/
1839 /*      */
1840 /*********************************COPYRIGHT NOTICE**********************************/