binaries/src/tcoffee/t_coffee_source/util_aln_analyze.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 #include "io_lib_header.h"
   8 #include "util_lib_header.h"
   9 #include "dp_lib_header.h"
  10 #include "define_header.h"
  11 /************************************************************************************/
  12 /*                NEW      ANALYZE 2    : SAR                                        */
  13 /************************************************************************************/
  14 float display_prediction_old (int **prediction, int n, Alignment *A, Alignment *S, int field);
  15
  16 float display_prediction (int ***count, Alignment *S, int c, int n);
  17 Alignment * filter_aln4sar0 ( Alignment *A, Alignment *S, int c, int leave, char *mode);
  18 Alignment * filter_aln4sar1 ( Alignment *A, Alignment *S, int c, int leave, char *mode);
  19 Alignment * filter_aln4sar2 ( Alignment *A, Alignment *S, int c, int leave, char *mode);
  20 Alignment * filter_aln4sar3 ( Alignment *A, Alignment *S, int c, int leave, char *mode);
  21 Alignment * filter_aln4sar4 ( Alignment *A, Alignment *S, int c, int leave, char *mode);
  22 Alignment * filter_aln4sar5 ( Alignment *A, Alignment *S, int c, int leave, char *mode);
  23
  24 int **sar2profile ( Alignment *A, Alignment *S, int c, int leave);
  25 int **sar2profile_sim ( Alignment *A, Alignment *S, int **sim, int comp, int leave);
  26 int sar_profile2score ( char *seq, int **profile);
  27 double sar_vs_iseq1( char *sar, int  *seq, float gl, int **sim, char *best_aa);
  28 double sar_vs_seq1 ( char *sar, char *seq, float gl, int **sim, char *best_aa);
  29 double sar_vs_seq2 ( char *sar, char *seq, float ng, int **mat, char *a);
  30 double sar_vs_seq3 ( char *sar, char *seq, float ng, int **mat, char *a);
  31
  32 double sar_vs_iseq4 ( char *sar, int *seq, float ng, int **mat, char *a);//supports an extended alphabet
  33 double sar_vs_seq4 ( char *sar, char *seq, float ng, int **mat, char *a);
  34
  35 double sar_vs_seq5 ( char *sar, char *seq, float ng, int **mat, char *a);
  36 int make_sim_pred ( Alignment *A,Alignment *S, int comp, int seq);
  37
  38 int **sar2profile_sim ( Alignment *A, Alignment *S, int **sim, int comp, int leave)
  39 {
  40
  41   int a, b, r, c, c1, c2, r1, r2, s, p;
  42   int ***cache, **profile;
  43
  44
  45   profile=declare_int (A->len_aln, 26);
  46   cache=declare_arrayN (3,sizeof (int),2,A->len_aln, 26);
  47
  48   for ( a=0; a< A->len_aln; a++)
  49     for ( b=0; b< A->nseq; b++)
  50       {
  51         r=tolower(A->seq_al[b][a]);
  52         c=( S->seq_al[comp][b]=='I')?1:0;
  53         if (b==leave || is_gap(r)) continue;
  54         cache [c][a][r-'a']++;
  55       }
  56   for (a=0; a< A->nseq; a++)
  57     {
  58       if ( a==leave) continue;
  59       for ( b=0; b< A->nseq; b++)
  60         {
  61           c1=(S->seq_al[comp][a]=='I')?1:0;
  62           c2=(S->seq_al[comp][b]=='I')?1:0;
  63           if ( b==leave || b==a || c1!=1 || c1==c2) continue;
  64           s=sim[a][b];
  65
  66           for (p=0; p<A->len_aln; p++)
  67             {
  68               r1=tolower(A->seq_al[a][p]);
  69               r2=tolower(A->seq_al[b][p]);
  70               if ( is_gap(r1) || is_gap(r2) || r1==r2)continue;
  71               r1-='a';r2-='a';
  72               if (cache[1][p][r2])continue;
  73               if ( s<50)continue;
  74               profile[p][r2]-=s;
  75             }
  76         }
  77     }
  78
  79   free_arrayN((void***)cache,3);
  80   return profile;
  81
  82 }
  83 int **sar2profile ( Alignment *A, Alignment *S, int comp, int leave)
  84 {
  85
  86   int a, b,c,r, n, v, npos=0;
  87   int ***cache, **profile;
  88   int ncat;
  89   float n_gap, max_gap;
  90   profile=declare_int (A->len_aln, 26);
  91   cache=declare_arrayN (3,sizeof (int),2,A->len_aln, 26);
  92
  93
  94
  95   for ( n=0, a=0; a< A->nseq; a++)
  96     {
  97       if ( a==leave) continue;
  98       else n+=(S->seq_al[comp][a]=='I')?1:0;
  99     }
 100
 101   for ( a=0; a< A->len_aln; a++)
 102     for ( b=0; b< A->nseq; b++)
 103       {
 104         r=tolower(A->seq_al[b][a]);
 105         c=( S->seq_al[comp][b]=='I')?1:0;
 106         if (b==leave) continue;
 107         else if (is_gap(r))continue;
 108         r-='a';
 109         cache [c][a][r]++;
 110       }
 111
 112   ncat=15; /*ncat: limit the analysis to columns containing less than ncat categories of aa*/
 113   max_gap=0.05;
 114   for (a=0; a< A->len_aln; a++)
 115     {
 116       for (n_gap=0,b=0; b< A->nseq; b++)
 117         n_gap+=(is_gap(A->seq_al[b][a]));
 118       n_gap/=(float)A->nseq;
 119
 120       if ( n_gap> max_gap)continue;
 121
 122       for (v=0,r=0; r< 26; r++)
 123         {
 124           if (cache [0][a][r] || cache[1][a][r])v++;
 125         }
 126
 127       for (n=0,r=0; r< 26 && v<ncat; r++)
 128         {
 129           if (cache [0][a][r] && !cache[1][a][r])
 130             {
 131               n++;
 132               profile[a][r]=-cache[0][a][r];
 133             }
 134         }
 135       if (n) npos++;
 136     }
 137
 138   free_arrayN((void***)cache,3);
 139   return profile;
 140
 141 }
 142 Alignment * filter_aln4sar0 ( Alignment *A, Alignment *S, int comp, int leave, char *mode)
 143 {
 144   return copy_aln (A,NULL);
 145 }
 146 Alignment * filter_aln4sar1 ( Alignment *inA, Alignment *S, int comp, int leave, char *mode)
 147 {
 148   Alignment *F, *A;
 149   int a, b,c, i,r, n0, n1,g, score;
 150   int ***cache, **list1, **list2;
 151   int Delta;
 152
 153   int T1;
 154
 155   /*Keep only the positions where there are residues ONLY associated with 0 sequences*/
 156
 157   list1=declare_int ( inA->nseq, 2);
 158   list2=declare_int ( inA->len_aln, 2);
 159
 160   cache=declare_arrayN (3,sizeof (int),inA->len_aln,2, 26);
 161   F=copy_aln (inA, NULL);
 162
 163   A=copy_aln (inA, NULL);
 164   A->nseq=strlen (S->seq_al[comp]);
 165
 166   strget_param (mode, "_T1_", "5", "%d", &T1);
 167   for ( a=0; a< A->len_aln; a++)
 168     {
 169       n1=n0=g=0;
 170       for (b=0; b< A->nseq; b++)
 171         {
 172           if ( b==leave) continue;
 173           i=(S->seq_al[comp][b]=='I')?1:0;
 174           r=tolower(A->seq_al[b][a]);
 175           if ( r=='-')continue;
 176           cache[a][i][r-'a']++;
 177         }
 178     }
 179
 180   for (a=0; a< A->nseq; a++)
 181     for ( score=0,b=0; b<A->len_aln; b++)
 182       {
 183         r=tolower (A->seq_al[a][b]);
 184         if ( is_gap(r))continue;
 185         else if ( cache[b][0][r-'a'] && !cache[b][1][r-'a'])list1[a][0]++;
 186       }
 187
 188   for (a=0; a< A->len_aln; a++)
 189     {
 190       for ( score=0,b=0; b< A->nseq; b++)
 191         {
 192           r=tolower (A->seq_al[b][a]);
 193           if ( r=='-')continue;
 194           else r-='a';
 195           if ( cache[a][0][r] && !cache[a][1][r])score ++;
 196         }
 197       list2[a][0]=a;
 198       list2[a][1]=score;
 199     }
 200   sort_int (list2, 2, 1, 0, F->len_aln-1);
 201
 202   Delta=A->len_aln/(100/T1);
 203   for ( a=0; a< F->len_aln-Delta; a++)
 204     {
 205       b=list2[a][0];
 206       for ( c=0; c<F->nseq; c++)
 207         {
 208           F->seq_al[c][b]='-';
 209         }
 210     }
 211
 212   ungap_aln (F);
 213   free_aln (A);
 214   free_arrayN ( (void ***)cache, 3);
 215   free_arrayN ((void**)list1, 2);
 216   free_arrayN ((void**)list2, 2);
 217
 218   return F;
 219 }
 220 Alignment * filter_aln4sar2 ( Alignment *inA, Alignment *S, int comp, int leave, char *mode)
 221 {
 222   Alignment *F, *A;
 223   int a,b,r,ncat;
 224   int *cache;
 225   int max_ncat=10;
 226
 227   /*Keep Low entropy columns that contain less than ncat categories of different amino acids*/
 228   /*REmove columns containing 10% or more gaps*/
 229
 230   cache=vcalloc ( 500, sizeof (char));
 231   F=copy_aln (inA, NULL);
 232   A=copy_aln (inA, NULL);
 233   A->nseq=strlen (S->seq_al[comp]);
 234   for ( a=0; a< A->len_aln; a++)
 235     {
 236       for (ncat=0,b=0; b< A->nseq; b++)
 237         {
 238           if ( b==leave) continue;
 239
 240           r=tolower(A->seq_al[b][a]);
 241           if ( !cache[r])ncat++;
 242           cache[r]++;
 243         }
 244
 245       if ( ncat <max_ncat && ((cache['-']*100)/A->nseq)<10)
 246         {
 247           ;
 248         }
 249       else
 250         {
 251           for (b=0; b<F->nseq; b++)
 252             {
 253               r=tolower(F->seq_al[b][a]);
 254               F->seq_al[b][a]='-';
 255               cache[r]=0;
 256             }
 257         }
 258       for (b=0; b<A->nseq; b++)
 259           {
 260             r=tolower(A->seq_al[b][a]);
 261             cache[r]=0;
 262           }
 263     }
 264
 265   free_aln (A);
 266   ungap_aln (F);
 267   vfree (cache);
 268   return F;
 269 }
 270
 271 Alignment * filter_aln4sar3 ( Alignment *inA, Alignment *S, int comp, int leave, char *mode)
 272 {
 273   Alignment *F, *rA, *A;
 274   int a, b,c;
 275   int **list1;
 276   char *bufS, *bufA;
 277   int Delta;
 278   int T3;
 279
 280   /*Keep the 10% positions most correlated with the 0/1 pattern*/
 281
 282   A=copy_aln (inA, NULL);
 283   A->nseq=strlen (S->seq_al[comp]);
 284   F=copy_aln (inA, NULL);
 285   rA=rotate_aln (A, NULL);
 286
 287   strget_param (mode, "_T3_", "10", "%d", &T3);
 288
 289
 290   list1=declare_int ( inA->len_aln, 2);
 291   bufA=vcalloc ( A->nseq+1, sizeof (char));
 292   bufS=vcalloc ( A->nseq+1, sizeof (char));
 293
 294   sprintf ( bufS, "%s", S->seq_al[comp]);
 295   splice_out_seg(bufS,leave, 1);
 296
 297
 298   for (a=0; a< A->len_aln; a++)
 299     {
 300       char aa;
 301       list1[a][0]=a;
 302       sprintf (bufA, "%s", rA->seq_al[a]);
 303       splice_out_seg (bufA,leave,1);
 304       list1[a][1]=(int)sar_vs_seq3 ( bufS, bufA,0,NULL, &aa);
 305     }
 306
 307   sort_int (list1, 2, 1, 0, F->len_aln-1);
 308   Delta=F->len_aln/(100/T3);
 309   for ( a=0; a< F->len_aln-Delta; a++)
 310     {
 311           b=list1[a][0];
 312
 313           for ( c=0; c<F->nseq; c++)
 314             {
 315               F->seq_al[c][b]='-';
 316             }
 317
 318     }
 319   F->score_aln=list1[F->len_aln-1][1];
 320   ungap_aln (F);
 321
 322   free_aln (rA);
 323   free_aln(A);
 324   free_arrayN ((void**)list1, 2);
 325   vfree (bufS);vfree (bufA);
 326   return F;
 327 }
 328 Alignment * filter_aln4sar4 ( Alignment *inA, Alignment *S, int comp, int leave, char *mode)
 329 {
 330   Alignment *F, *A;
 331   int a, b,c, i,r, n0, n1,g,score;
 332   int ***cache, **list1, **list2;
 333
 334   /*Keep only the positions where there are residues ONLY associated with 0 sequences*/
 335
 336   list1=declare_int ( inA->nseq, 2);
 337   list2=declare_int ( inA->len_aln, 2);
 338
 339   cache=declare_arrayN (3,sizeof (int),inA->len_aln,2, 26);
 340   F=copy_aln (inA, NULL);
 341   A=copy_aln (inA, NULL);
 342   A->nseq=strlen (S->seq_al[comp]);
 343
 344   for ( a=0; a< A->len_aln; a++)
 345     {
 346       n1=n0=g=0;
 347       for (b=0; b< A->nseq; b++)
 348         {
 349           if ( b==leave) continue;
 350           i=(S->seq_al[comp][b]=='I')?1:0;
 351           r=tolower(A->seq_al[b][a]);
 352           if ( r=='-')continue;
 353           cache[a][i][r-'a']++;
 354           n1+=i;
 355         }
 356     }
 357
 358
 359   for (a=0; a< A->len_aln; a++)
 360     {
 361       for ( score=0,b=0; b< A->nseq; b++)
 362         {
 363           r=tolower (F->seq_al[b][a]);
 364           if ( r=='-')continue;
 365           else r-='a';
 366           if (cache[a][1][r]>=n1/2)score=1;
 367         }
 368       list2[a][0]=a;
 369       list2[a][1]=score;
 370     }
 371
 372
 373   for ( a=0; a< F->len_aln; a++)
 374     {
 375       if ( list2[a][1]==1);
 376       else
 377         {
 378           b=list2[a][0];
 379           for ( c=0; c<F->nseq; c++)
 380             {
 381               F->seq_al[c][b]='-';
 382             }
 383         }
 384     }
 385   ungap_aln (F);
 386   free_aln (A);
 387   free_arrayN ( (void ***)cache, 3);
 388   free_arrayN ((void**)list1, 2);
 389   free_arrayN ((void**)list2, 2);
 390
 391   return F;
 392 }
 393
 394 Alignment * filter_aln4sar5 ( Alignment *inA, Alignment *S, int comp, int leave, char *mode)
 395 {
 396   Alignment *F, *rA, *A;
 397   int a, b,c;
 398   int **list1;
 399   char *bufS, *bufA;
 400   int max;
 401   /*Look for the positions that show the best correlation between the sequence variation and the SAR*/
 402
 403   A=copy_aln (inA, NULL);
 404   A->nseq=strlen (S->seq_al[comp]);
 405
 406   rA=rotate_aln (inA, NULL);
 407   F=copy_aln (inA, NULL);
 408
 409   list1=declare_int ( A->len_aln, 2);
 410   bufA=vcalloc ( A->nseq+1, sizeof (char));
 411   bufS=vcalloc ( A->nseq+1, sizeof (char));
 412
 413
 414
 415   sprintf ( bufS, "%s", S->seq_al[comp]);
 416   splice_out_seg(bufS,leave, 1);
 417
 418
 419   for (a=0; a< A->len_aln; a++)
 420     {
 421       char aa;
 422       list1[a][0]=a;
 423       sprintf (bufA, "%s", rA->seq_al[a]);
 424       splice_out_seg (bufA,leave,1);
 425       list1[a][1]=(int)sar_vs_seq4 ( bufS, bufA,0,NULL, &aa);
 426     }
 427
 428   sort_int (list1, 2, 1, 0, F->len_aln-1);
 429   max=F->score=list1[F->len_aln-1][1];
 430   max-=(max/10);
 431
 432
 433   for ( a=0; a< F->len_aln-10; a++)
 434     {
 435
 436           b=list1[a][0];
 437
 438           for ( c=0; c<F->nseq; c++)
 439             {
 440               F->seq_al[c][b]='-';
 441             }
 442
 443     }
 444   F->score_aln=10;
 445   ungap_aln (F);
 446   free_aln (inA);
 447   free_aln (rA);
 448   free_arrayN ((void**)list1, 2);
 449   vfree (bufS);vfree (bufA);
 450   return F;
 451 }
 452
 453 int sar_profile2score ( char *seq, int **P)
 454 {
 455   int a,r, l, score;
 456
 457   l=strlen (seq);
 458   for ( score=0,a=0; a< l; a++)
 459     {
 460       r=seq[a];
 461       if ( is_gap(r))continue;
 462       score+=P[a][tolower(r)-'a'];
 463     }
 464   return score;
 465 }
 466 int make_sim_pred ( Alignment *A,Alignment *S, int comp, int seq)
 467 {
 468   int a, b, i, r1, r2;
 469   static float **cscore;
 470   static float **tscore;
 471
 472   if ( !cscore)
 473     {
 474       cscore=declare_float (2, 2);
 475       tscore=declare_float (2, 2);
 476     }
 477
 478   for (a=0; a< 2; a++)for (b=0; b<2; b++)cscore[a][b]=tscore[a][b]=0;
 479
 480   for ( a=0; a<A->len_aln; a++)
 481     {
 482       r1=A->seq_al[seq][a];
 483       if ( r1=='-') continue;
 484       else
 485         {
 486           for ( b=0; b< A->nseq; b++)
 487             {
 488               if (b==seq) continue;
 489               else
 490                 {
 491                   r2=A->seq_al[b][a];
 492                   if (r2=='-')continue;
 493                   else
 494                     {
 495
 496                       i=(S->seq_al[comp][b]=='I')?1:0;
 497                       cscore[i][0]+=(r1==r2)?1:0;
 498                       cscore[i][1]++;
 499                     }
 500                 }
 501             }
 502
 503           for (i=0; i<2; i++)
 504             {
 505               cscore[i][0]/=(cscore[i][1]==0)?1:cscore[i][1];
 506               tscore[i][0]+=cscore[i][0];tscore[i][1]++;
 507               cscore[i][0]=cscore[i][1]=0;
 508             }
 509         }
 510     }
 511
 512   fprintf ( stdout, "\nn\t 1: %.2f 0: %.2f", tscore[1][0],tscore[0][0]);
 513   return ( tscore[1][0]>=tscore[0][0])?1:0;
 514 }
 515
 516
 517 Alignment * sar_analyze (Alignment *inA, Alignment *inS, char *mode)
 518 {
 519   int ***sim,***glob_results, ***comp_results;
 520   int *count;
 521   int a,b,c,m;
 522   float *tot2;
 523   Alignment *A=NULL,*S=NULL,*F, *SUBSET;
 524   char *subset, *target;
 525   int jack, T, filter;
 526   filter_func *ff;
 527   int n_methods=0;
 528   char *prediction, *reliability;
 529   int pred_start=0, pred_end, ref_start=0, ref_end;
 530   int display, CSV=1, NONCSV=0;
 531   char method[5];
 532
 533   strget_param (mode, "_METHOD_", "1111", "%s_", method);
 534   ff=vcalloc (6,sizeof (filter_func));
 535   if (method[0]=='1')ff[n_methods++]=filter_aln4sar0;
 536   if (method[1]=='1')ff[n_methods++]=filter_aln4sar1;
 537   if (method[2]=='1')ff[n_methods++]=filter_aln4sar2;
 538   if (method[3]=='1')ff[n_methods++]=filter_aln4sar3;
 539   /*
 540     ff[n_methods++]=filter_aln4sar4;
 541     ff[n_methods++]=filter_aln4sar5;
 542   */
 543   sim=vcalloc (n_methods, sizeof (int**));
 544
 545
 546   tot2=vcalloc ( 10, sizeof (float));
 547   subset=vcalloc ( 100, sizeof (char));
 548   target=vcalloc ( 100, sizeof (char));
 549
 550   strget_param (mode, "_TARGET_", "no", "%s_", target);
 551   strget_param (mode, "_SUBSET_", "no", "%s_", subset);
 552   strget_param (mode, "_JACK_", "0", "%d", &jack);
 553   strget_param (mode, "_T_", "0", "%d", &T);
 554   strget_param (mode, "_FILTER_", "11", "%d", &filter);
 555   strget_param (mode, "_DISPLAY_", "0", "%d", &display);
 556
 557
 558
 559   if ( !strm (target, "no"))
 560     {
 561       Alignment *T;
 562       T=main_read_aln(target, NULL);
 563       if ( T->len_aln !=inA->len_aln )
 564         {
 565           printf_exit ( EXIT_FAILURE,stderr, "Error: %s is incompatible with the reference alignment [FATAL:%s]",target,PROGRAM);
 566         }
 567
 568       inA=stack_aln (inA, T);
 569
 570     }
 571
 572   if ( !strm(subset, "no"))
 573     {
 574       SUBSET=main_read_aln (subset, NULL);
 575       sarset2subsarset ( inA, inS, &A, &S, SUBSET);
 576     }
 577   else
 578     {
 579       A=inA;
 580       S=inS;
 581     }
 582
 583
 584   prediction=vcalloc ( n_methods+1, sizeof (char));
 585   reliability=vcalloc ( n_methods+1, sizeof (char));
 586
 587   glob_results=declare_arrayN(3, sizeof (int), n_methods*2, 2, 2);
 588
 589   count=vcalloc (S->nseq, sizeof (int));
 590   for (a=0; a<S->nseq; a++)
 591     {
 592       int l;
 593       l=strlen (S->seq_al[a]);
 594       for ( b=0; b<l; b++)
 595         count[a]+=(S->seq_al[a][b]=='I')?1:0;
 596     }
 597   if ( display==CSV)
 598     {fprintf ( stdout, "\nCompound %s ; Ntargets %d", S->name[a],count[a]);
 599       pred_start=(strlen (S->seq_al[0])==A->nseq)?0:strlen (S->seq_al[0]);
 600       pred_end=A->nseq;
 601       for (a=pred_start; a< pred_end; a++)
 602         fprintf ( stdout, ";%s", A->name[a]);
 603       fprintf ( stdout, ";npred;");
 604     }
 605
 606
 607   for (a=0; a<S->nseq; a++)
 608     {
 609       int n_pred;
 610       comp_results=declare_arrayN(3, sizeof (int), n_methods*2, 2, 2);
 611
 612       pred_start=(strlen (S->seq_al[a])==A->nseq)?0:strlen (S->seq_al[a]);
 613       pred_end=A->nseq;
 614       if ( display==CSV)fprintf ( stdout, "\n%s;%d", S->name[a],count[a]);
 615
 616       for (n_pred=0,b=pred_start; b<pred_end;b++)
 617         {
 618           int t, score=0,pred, real;
 619
 620           if ( display==NONCSV)fprintf ( stdout, "\n>%-15s %10s %c ", S->name[a], A->name[b], (pred_start==0)?S->seq_al[a][b]:'?');
 621           if (jack || b==pred_start)
 622             {
 623               for (m=0; m<n_methods; m++)
 624                 {
 625                   free_int (sim[m], -1);
 626                   F=(ff[m]) (A,S,a,(jack==0)?-1:b, mode);
 627                   sim[m]=aln2sim_mat(F, "idmat");
 628                   free_aln (F);
 629                 }
 630             }
 631
 632           for (m=0; m<n_methods; m++)
 633             {
 634               int Nbsim=0,Ybsim=0,bsim=0;
 635               ref_start=0;
 636               ref_end=strlen (S->seq_al[m]);
 637
 638               for (c=ref_start;c<ref_end; c++)
 639                 {
 640                   if ( b==c) continue;
 641                   else if ( S->seq_al[a][c]=='O')
 642                     {
 643                       Nbsim=MAX(Nbsim,sim[m][b][c]);
 644                     }
 645                   else
 646                     {
 647                       Ybsim=MAX(Ybsim,sim[m][b][c]);
 648                     }
 649                 }
 650
 651               bsim=(Ybsim>Nbsim)?Ybsim:-Nbsim;
 652               pred=(bsim>0)?1:0;
 653               real=(S->seq_al[a][b]=='O')?0:1;
 654               comp_results[m][pred][real]++;
 655               glob_results[m][pred][real]++;
 656               score+=pred;
 657               prediction[m]=pred+'0';
 658               reliability[m]=(FABS((Ybsim-Nbsim))-1)/10+'0';
 659             }
 660
 661           if ( score>0)n_pred++;
 662           prediction[m]=reliability[m]='\0';
 663           if (display==NONCSV)fprintf ( stdout, "Compound_Count:%d primary_predictions: %s Total: %d", count[a],prediction, score);
 664           else if ( display==CSV)fprintf ( stdout, ";%d", score);
 665           for (t=0; t<n_methods; t++)
 666             {
 667               if (score>t)
 668                 {
 669                   comp_results[t+n_methods][1][real]++;
 670                   glob_results[t+n_methods][1][real]++;
 671                 }
 672               else
 673                 {
 674                   comp_results[t+n_methods][0][real]++;
 675                   glob_results[t+n_methods][0][real]++;
 676                 }
 677             }
 678         }
 679       if ( display==NONCSV)
 680         {if ( pred_start==0)display_prediction (comp_results, S,a, n_methods*2);}
 681       else fprintf (stdout, ";%d;",n_pred);
 682     }
 683   if ( display==NONCSV)if (pred_start==0)display_prediction (glob_results, S,-1, n_methods*2);
 684
 685
 686   exit (EXIT_SUCCESS);
 687 }
 688 float display_prediction (int ***count, Alignment *S, int c, int n)
 689 {
 690   float tp,tn,fn,fp,sp,sn,sn2;
 691   int a, nm;
 692
 693   nm=n/2;
 694
 695   for (a=0; a<n; a++)
 696     {
 697       tp=count[a][1][1];
 698       tn=count[a][0][0];
 699       fp=count[a][1][0];
 700       fn=count[a][0][1];
 701
 702       sn2=tp/(tp+fp);
 703       sn=tp/(tp+fn);
 704       sp=tn/(tn+fp);
 705       if ( a<nm)fprintf ( stdout, "\n>#Method %d Compound %15s sp=%.2f sn=%.2f sn2=%.2f",a, (c==-1)?"TOTAL":S->name[c],sp, sn, sn2 );
 706       else fprintf ( stdout, "\n>#Combined: T=%d Compound %15s sp=%.2f sn=%.2f sn2=%.2f",a-nm, (c==-1)?"TOTAL":S->name[c],sp, sn, sn2 );
 707     }
 708   fprintf ( stdout, "\n");
 709   return 0;
 710 }
 711
 712 float display_prediction_2 (int **prediction, int n,Alignment *A, Alignment *S, int field)
 713 {
 714   int a, t, T;
 715   float max_sn, max_sp;
 716
 717   if ( field==17 || field ==18)
 718     {
 719       printf_exit ( EXIT_FAILURE, stderr, "\nERROR: Do not use filed %d in display_prediction", field);
 720     }
 721
 722   sort_int_inv ( prediction, 10,field, 0, n-1);
 723   for (t=0,a=0; a<n; a++)
 724     {
 725       t+=prediction[a][3];
 726       prediction[a][17]=t;
 727     }
 728
 729   for (t=0,a=n-1; a>=0; a--)
 730     {
 731       prediction[a][18]=t;
 732       t+=prediction[a][3];
 733     }
 734
 735   max_sn=max_sp=T=0;
 736   for (a=0; a<n; a++)
 737     {
 738       float tp, fn, fp, sp, sn;
 739
 740       tp=prediction[a][17];
 741       fn=prediction[a][18];
 742       fp=(a+1)-tp;
 743
 744       sp=((tp+fp)==0)?0:tp/(tp+fp);
 745       sn=((tp+fn)==0)?0:tp/(tp+fn);
 746
 747       if (sp>0.8)
 748         {
 749           if (sn>max_sn)
 750             {
 751               max_sn=sn;
 752               max_sp=sp;
 753
 754               T=prediction[a][field];
 755             }
 756         }
 757     }
 758   if (max_sn==0)
 759       fprintf (stdout, "\n T =%d SN=%.2f SP= %.2f",T,max_sn,max_sp);
 760   else
 761       fprintf (stdout, "\n T =%d SN=%.2f SP= %.2f",T,max_sn,max_sp);
 762
 763   return max_sn;
 764 }
 765
 766
 767 /************************************************************************************/
 768 /*                NEW      ANALYZE     : SAR                                        */
 769 /************************************************************************************/
 770 float** cache2pred1 (Alignment *A,int**cache, int *ns, int **ls, Alignment *S, char *compound, char *mode);
 771 float** cache2pred2 (Alignment *A,int**cache, int *ns, int **ls, Alignment *S, char *compound, char *mode);
 772 float** cache2pred3 (Alignment *A,int**cache, int *ns, int **ls, Alignment *S, char *compound, char *mode);
 773 float** cache2pred4 (Alignment *A,int**cache, int *ns, int **ls, Alignment *S, char *compound, char *mode);
 774 float** cache2pred5 (Alignment *A,int**cache, int *ns, int **ls, Alignment *S, char *compound, char *mode);
 775 float** cache2pred_new (Alignment *A,int**cache, int *ns, int **ls, Alignment *S, char *compound, char *mode);
 776
 777 int **sar2cache_adriana ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode);
 778 int **sar2cache_proba_old ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode);
 779 int **sar2cache_count1 ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode);
 780 int **sar2cache_count2 ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode);
 781 int **sar2cache_count3 ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode);
 782
 783 int **sar2cache_proba_new ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode);
 784 int **sar2cache_proba2 ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode);
 785 int **analyze_sar_compound1 ( char *name, char *seq, Alignment *A, char *mode);
 786 int **analyze_sar_compound2 ( char *name, char *seq, Alignment *A, char *mode);
 787
 788 int aln2n_comp_col ( Alignment *A, Alignment *S, int ci);
 789
 790
 791
 792
 793 int ***simple_sar_analyze_vot ( Alignment *inA, Alignment *SAR, char *mode);
 794 int ***simple_sar_analyze_col ( Alignment *inA, Alignment *SAR, char *mode);
 795
 796
 797
 798 int sarset2subsarset ( Alignment *A, Alignment *S, Alignment **subA, Alignment **subS, Alignment *SUB);
 799 int benchmark_sar (int v);
 800 int aln2jack_group1 (Alignment *A, int seq, int **l1, int *nl1, int **l2, int *nl2);
 801 int aln2jack_group2 (Alignment *A, int seq, int **l1, int *nl1, int **l2, int *nl2);
 802 int aln2jack_group3 (Alignment *A, char *sar_seq, int **l1, int *nl1, int **l2, int *nl2);
 803 float** jacknife5 (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode);
 804 float** jacknife6 (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode);
 805
 806 int process_cache ( Alignment *A,Alignment *S, int ***Cache, char *mode);
 807 Alignment *analyze_compounds (Alignment *A, Alignment *S, char *mode);
 808
 809 Alignment *analyze_compounds (Alignment *A, Alignment *S, char *mode)
 810 {
 811   int a, b, c, tot, n;
 812   int **sim;
 813   int sar1, sar2;
 814
 815   sim=aln2sim_mat (A, "idmat");
 816   for (a=0; a< S->nseq; a++)
 817     {
 818       for (n=0, tot=0, b=0; b< A->nseq-1; b++)
 819         {
 820           sar1=(S->seq_al[a][b]=='I')?1:0;
 821           for ( c=b+1; c<A->nseq; c++)
 822             {
 823               sar2=(S->seq_al[a][c]=='I')?1:0;
 824
 825               if (sar1 && sar2)
 826                 {
 827                   tot+=sim[b][c];
 828                   n++;
 829                 }
 830             }
 831         }
 832       fprintf ( stdout, ">%-10s   CMPSIM: %.2f\n", S->name[a],(float)tot/(float)n);
 833     }
 834   free_int (sim, -1);
 835   return A;
 836 }
 837
 838 int print_seq_pos ( int pos, Alignment *A, char *seq);
 839 int abl1_evaluation (int p);
 840 int print_seq_pos ( int pos, Alignment *A, char *seq)
 841 {
 842   int a, b, s;
 843
 844   s=name_is_in_list (seq, A->name, A->nseq, MAXNAMES);
 845   fprintf ( stdout, "S=%d", s);
 846
 847   for (b=0,a=0; a<pos; a++)
 848     {
 849       if (!is_gap (A->seq_al[s][a]))b++;
 850     }
 851   fprintf ( stdout, "Pos %d SEQ %s: %d ", pos+1, seq, b+246);
 852   if ( strm ( seq, "ABL1")) fprintf ( stdout , "PT: %d", abl1_evaluation (b+246));
 853   return 0;
 854 }
 855
 856 int process_cache ( Alignment *A,Alignment *S, int  ***Cache, char *mode)
 857 {
 858   int a, b;
 859   int **pos, **pos2;
 860   int **C;
 861   int ab1, *ab1_pos;
 862   int weight_mode;
 863
 864   strget_param ( mode, "_WEIGHT_", "1", "%d", &weight_mode);
 865   pos=declare_int(A->len_aln+1,2);
 866   pos2=declare_int (A->len_aln+1,S->nseq);
 867   for (a=0; a<S->nseq; a++)
 868     {
 869       C=Cache[a];
 870       for (b=0; b< A->len_aln; b++)
 871         {
 872             pos[b][0]+=C[26][b];
 873             if ( C[26][b]>0)
 874               {
 875                 pos[b][1]++;
 876                 pos2[b][a]=1;
 877               }
 878         }
 879     }
 880
 881   C=Cache[0];
 882   ab1=name_is_in_list ("ABL1", A->name, A->nseq,100);
 883   ab1_pos=vcalloc (A->len_aln+1, sizeof (int));
 884
 885   for ( b=0,a=0; a< A->len_aln; a++)
 886     {
 887       if ( A->seq_al[ab1][a]=='-')ab1_pos[a]=-1;
 888       else ab1_pos[a]=++b;
 889     }
 890
 891   for ( a=0; a< A->len_aln; a++)
 892     {
 893       fprintf ( stdout, "\n%4d %5d %5d %5d [%c] [%2d] ALN", a+1, pos[a][0], pos[a][1], ab1_pos[a]+246,A->seq_al[ab1][a],abl1_evaluation (ab1_pos[a]+246));
 894       for ( b=0; b< S->nseq; b++)fprintf ( stdout, "%d", pos2[a][b]);
 895     }
 896   return 1;
 897 }
 898 int abl1_evaluation (int p)
 899 {
 900   if ( p==248) return 10;
 901   if ( p==250) return 10;
 902   if ( p==253) return 10;
 903   if ( p==254) return 10;
 904   if ( p==255) return 9;
 905   if ( p==256) return 10;
 906   if ( p==257) return 5;
 907   if ( p==258) return 8;
 908   if ( p==269) return 8;
 909   if ( p==291) return 4;
 910   if ( p==294) return 8;
 911   if ( p==299) return 10;
 912   if ( p==306) return 0;
 913   if ( p==314) return 9;
 914   if ( p==315) return 10;
 915   if ( p==318) return 10;
 916
 917   if ( p==319) return 10;
 918   if ( p==321) return 10;
 919   if ( p==323) return 0;
 920   if ( p==324) return 0;
 921   if ( p==339) return 0;
 922   if ( p==340) return 0;
 923   if ( p==355) return 5;
 924   if ( p==364) return 10;
 925
 926   if ( p==366) return 0;
 927   if ( p==368) return 10;
 928   if ( p==370) return 10;
 929   if ( p==372) return 0;
 930   if ( p==378) return 8;
 931   if ( p==382) return 10;
 932
 933   if ( p==384) return 10;
 934   if ( p==387) return 10;
 935   if ( p==395) return 8;
 936
 937   if ( p==398) return 8;
 938   if ( p==399) return 8;
 939   if ( p==400) return 8;
 940   if ( p==403) return 0;
 941   if ( p==416) return 8;
 942   if ( p==419) return 5;
 943   if ( p>400) return 0;
 944   return -1;
 945 }
 946 float** cache2pred1 (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode)
 947 {
 948   int s1, s2, seq1, seq2, r1, r2,col, pred, real, ci;
 949   double score, max, id, m;
 950   float **R, T;
 951
 952
 953   int used_col, used_res,is_used_col, n_res=0;
 954   int weight_mode;
 955   /*Predict on ns[1] what was trained on ns[0]*/
 956
 957   strget_param ( mode, "_THR_", "0.09", "%f", &T);
 958   strget_param ( mode, "_WEIGHT_", "0", "%d", &weight_mode);
 959
 960   R=declare_float (2, 2);
 961   ci=name_is_in_list ( compound, S->name, S->nseq, -1);
 962
 963
 964
 965   for (s1=0; s1<ns[1]; s1++)
 966     {
 967       int v;
 968       seq1=ls[1][s1];
 969
 970       for (max=0,score=0, col=0; col<A->len_aln; col++)
 971         {
 972           int max1;
 973           r1=tolower (A->seq_al[seq1][col]);
 974           for (max1=0,id=0, m=0,s2=0; s2<ns[0]; s2++)
 975             {
 976               seq2=ls[0][s2];
 977               if ( S->seq_al[ci][seq2]=='O')continue;
 978               if ( cache[seq2][col]==0 && !is_gap( A->seq_al[seq2][col]))continue;
 979
 980               r2=tolower ( A->seq_al[seq2][col]);
 981               if ( is_gap(r2))continue;
 982
 983               v=(cache[seq2][col]>0 && weight_mode==1)?cache[seq2][col]:1;
 984
 985               max+=v;
 986               if ( r2==r1)
 987                 {
 988                   score+=v;
 989                 }
 990
 991             }
 992
 993         }
 994       pred=(( score/max) >T)?1:0;
 995       real=(S->seq_al[ci][seq1]=='I')?1:0;
 996       R[pred][real]++;
 997
 998       fprintf ( stdout, "\n>%s %d%d SCORE %.2f C %s [SEQ]\n", A->name[seq1],real, pred, (float)score/(float)max, compound);
 999     }
1000
1001   for (used_col=0,used_res=0,col=0; col<A->len_aln; col++)
1002     {
1003       for (is_used_col=0,s2=0; s2<ns[0]; s2++)
1004         {
1005           seq2=ls[0][s2];
1006           if ( cache[seq2][col]==0 && !is_gap(A->seq_al[seq2][col]))n_res++;
1007           else if (is_gap(A->seq_al[seq2][col]));
1008           else
1009             {
1010             is_used_col=1;
1011             used_res++;
1012             }
1013         }
1014       used_col+=is_used_col;
1015     }
1016   fprintf ( stdout, "\n>%s USED_POSITIONS: COL: %.2f RES: %.2f COMP\n", S->name[ci],  (float)used_col/(float)A->len_aln, (float)used_res/(float) n_res);
1017
1018   return R;
1019 }
1020
1021 float** cache2pred2 (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode)
1022 {
1023   int s1, s2, seq1, seq2, r1, r2,col, pred, real, ci;
1024   double score, max;
1025   float **R, T;
1026
1027
1028   int used_col, used_res,is_used_col, n_res=0;
1029   /*Predict on ns[1] what was trained on ns[0]*/
1030
1031   strget_param ( mode, "_THR_", "0.5", "%f", &T);
1032
1033
1034   R=declare_float (2, 2);
1035   ci=name_is_in_list ( compound, S->name, S->nseq, -1);
1036
1037   for (s1=0; s1<ns[1]; s1++)
1038     {
1039       int v;
1040       seq1=ls[1][s1];
1041       fprintf ( stdout, "\n");
1042       for (max=0,score=0, col=0; col<A->len_aln; col++)
1043         {
1044           int used;
1045
1046           r1=tolower (A->seq_al[seq1][col]);
1047           for (used=0,s2=0; s2<ns[0]; s2++)
1048             {
1049               seq2=ls[0][s2];
1050
1051               if ( S->seq_al[ci][seq2]=='O')continue;
1052               if ( cache[seq2][col]==0 && !is_gap( A->seq_al[seq2][col]))continue;
1053
1054
1055               r2=tolower ( A->seq_al[seq2][col]);
1056               if ( is_gap(r2))continue;
1057
1058               v=cache[seq2][col];
1059               if ( r2==r1){score+=v;}
1060               used=1;
1061               max+=v;
1062             }
1063           if (used) fprintf ( stdout, "%c", r1);
1064         }
1065
1066       pred=(( score/max) >T)?1:0;
1067       real=(S->seq_al[ci][seq1]=='I')?1:0;
1068       R[pred][real]++;
1069
1070       fprintf ( stdout, "PSEQ: %-10s SC: %4d MAX: %4d S: %.2f R: %4d", A->name[seq1],(int)score, (int)max, (float)score/max,real);
1071
1072     }
1073
1074   for (used_col=0,used_res=0,col=0; col<A->len_aln; col++)
1075     {
1076       for (is_used_col=0,s2=0; s2<ns[0]; s2++)
1077         {
1078           seq2=ls[0][s2];
1079           if ( cache[seq2][col]==0 && !is_gap(A->seq_al[seq2][col]))n_res++;
1080           else if (is_gap(A->seq_al[seq2][col]));
1081           else
1082             {
1083             is_used_col=1;
1084             used_res++;
1085             }
1086         }
1087       used_col+=is_used_col;
1088     }
1089   fprintf ( stdout, "\n>%s USED_POSITIONS: COL: %.2f RES: %.2f COMP\n", S->name[ci],  (float)used_col/(float)A->len_aln, (float)used_res/(float) n_res);
1090
1091   return R;
1092 }
1093
1094 float** cache2pred3 (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode)
1095 {
1096   int s1, s2, seq1, seq2, r1, r2,col, pred, real, ci, a, n;
1097   double score, max;
1098   float **R, T;
1099
1100
1101
1102   int tp, tn, fn, fp;
1103   int best_tp, best_fp;
1104   int delta, best_delta;
1105   int **list;
1106
1107   /*Predict on ns[1] what was trained on ns[0]*/
1108
1109   strget_param ( mode, "_THR_", "0.5", "%f", &T);
1110
1111
1112   R=declare_float (2, 2);
1113   ci=name_is_in_list ( compound, S->name, S->nseq, -1);
1114   list=declare_int ( ns[1],3);
1115
1116   for (s1=0; s1<ns[1]; s1++)
1117     {
1118       int v;
1119       seq1=ls[1][s1];
1120
1121       for (max=0,score=0, col=0; col<A->len_aln; col++)
1122         {
1123           int used;
1124
1125           r1=tolower (A->seq_al[seq1][col]);
1126           for (used=0,s2=0; s2<ns[0]; s2++)
1127             {
1128               seq2=ls[0][s2];
1129
1130               if ( S->seq_al[ci][seq2]=='O')continue;
1131               if ( cache[seq2][col]==0 && !is_gap( A->seq_al[seq2][col]))continue;
1132
1133
1134               r2=tolower ( A->seq_al[seq2][col]);
1135               if ( is_gap(r2))continue;
1136
1137               v=cache[seq2][col];
1138               if ( r2==r1){score+=v;}
1139               used=1;
1140               max+=v;
1141             }
1142         }
1143
1144
1145
1146       pred=(( score/max) >T)?1:0;
1147       real=(S->seq_al[ci][seq1]=='I')?1:0;
1148
1149       list[s1][0]=real;
1150       list[s1][1]=(int)((score/max)*(float)1000);
1151       list[s1][2]=seq1;
1152
1153
1154
1155     }
1156   sort_int_inv (list, 3, 1, 0, ns[1]-1);
1157
1158   for ( a=0; a<ns[1]; a++)
1159     {
1160       seq1=list[a][2];
1161       fprintf ( stdout, "PSEQ: %-10s SC: %5d R: %4d\n", A->name[seq1],list[a][0], list[a][1]);
1162     }
1163
1164   for (n=0, a=0; a<ns[1]; a++)n+=list[a][0];
1165   for (best_delta=100000,best_tp=0,tp=0,fp=0,best_fp=0,a=0; a<ns[1]; a++)
1166     {
1167       tp+=list[a][0];
1168       fp+=1-list[a][0];
1169       delta=(n-(tp+fp));
1170       if (FABS(delta)<best_delta)
1171         {
1172           best_delta=delta;
1173           best_tp=tp;
1174           best_fp=fp;
1175         }
1176     }
1177   /*R[pred][real]*/
1178   tp=best_tp;
1179   fp=best_fp;
1180   fn=n-tp;
1181   tn=ns[1]-(tp+fp+fn);
1182   R[1][1]=tp;
1183   R[1][0]=fp;
1184   R[0][1]=fn;
1185   R[0][0]=tn;
1186   free_int (list, -1);
1187   return R;
1188 }
1189 float** cache2pred4 (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode)
1190 {
1191   int s1, s2, seq1, seq2, ci, a,b, c, n;
1192   double score;
1193   float **R;
1194
1195
1196   int tp, tn, fn, fp;
1197   int best_tp, best_fp;
1198   int delta, best_delta;
1199   int **list;
1200   int **sim;
1201   int *ul;
1202   int nused=0;
1203
1204   /*Predict on ns[1] what was trained on ns[0]*/
1205   /*Identify interesting coloumns*/
1206   ul=vcalloc ( A->len_aln, sizeof (int));
1207   for (a=0; a< A->len_aln; a++)
1208     for ( b=0; b< A->nseq; b++)
1209       if ( cache[b][a])ul[nused++]=a;
1210
1211   /*compute the similarity on the used columns*/
1212
1213   R=declare_float (2, 2);
1214   sim=declare_int (A->nseq, A->nseq);
1215   for (a=0; a< A->nseq; a++)
1216     for ( b=0; b< A->nseq; b++)
1217       {
1218         for (c=0; c< nused; c++)
1219           {
1220             if ( A->seq_al[a][ul[c]]==A->seq_al[b][ul[c]])sim[a][b]++;
1221           }
1222         sim[a][b]=(sim[a][b]*100)/nused;
1223       }
1224   vfree (ul);
1225
1226
1227
1228
1229   ci=name_is_in_list ( compound, S->name, S->nseq, -1);
1230   list=declare_int ( ns[1],2);
1231
1232   for (s1=0; s1<ns[1]; s1++)
1233     {
1234
1235       seq1=ls[1][s1];
1236
1237       for (score=0,s2=0; s2<ns[0]; s2++)
1238         {
1239           seq2=ls[0][s2];
1240
1241           if ( seq1==seq2)continue;
1242           if (S->seq_al[ci][seq2]=='I')score=MAX(score, sim[seq1][seq2]);
1243         }
1244       list[s1][0]=(S->seq_al[ci][seq1]=='I')?1:0;
1245       list[s1][1]=(int)score;
1246
1247     }
1248   sort_int_inv (list, 2, 1, 0, ns[1]-1);
1249
1250   for (n=0, a=0; a<ns[1]; a++)n+=list[a][0];
1251   for (best_delta=100000,best_tp=0,tp=0,fp=0,best_fp=0,a=0; a<ns[1]; a++)
1252     {
1253       tp+=list[a][0];
1254       fp+=1-list[a][0];
1255       delta=(n-(tp+fp));
1256       if (FABS(delta)<best_delta)
1257         {
1258           best_delta=delta;
1259           best_tp=tp;
1260           best_fp=fp;
1261         }
1262     }
1263   /*R[pred][real]*/
1264   tp=best_tp;
1265   fp=best_fp;
1266   fn=n-tp;
1267   tn=ns[1]-(tp+fp+fn);
1268   R[1][1]=tp;
1269   R[1][0]=fp;
1270   R[0][1]=fn;
1271   R[0][0]=tn;
1272   free_int (list, -1);
1273   free_int (sim, -1);
1274   return R;
1275 }
1276
1277 float** cache2pred5 (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode)
1278 {
1279   int s1, s2, seq1, seq2, ci, a, n;
1280   double score;
1281   float **R;
1282
1283
1284
1285   int tp, tn, fn, fp;
1286   int best_tp, best_fp;
1287   int delta, best_delta;
1288   int **list;
1289   static int **sim;
1290
1291   /*Predict on ns[1] what was trained on ns[0]*/
1292
1293   R=declare_float (2, 2);
1294
1295   if ( sim==NULL)
1296     sim=aln2sim_mat (A, "idmat");
1297
1298
1299
1300   ci=name_is_in_list ( compound, S->name, S->nseq, -1);
1301   list=declare_int ( ns[1],2);
1302
1303   for (s1=0; s1<ns[1]; s1++)
1304     {
1305
1306       seq1=ls[1][s1];
1307
1308       for (score=0,s2=0; s2<ns[0]; s2++)
1309         {
1310           seq2=ls[0][s2];
1311
1312           if ( seq1==seq2)continue;
1313           if (S->seq_al[ci][seq2]=='I')score=MAX(score, sim[seq1][seq2]);
1314         }
1315       list[s1][0]=(S->seq_al[ci][seq1]=='I')?1:0;
1316       list[s1][1]=(int)score;
1317
1318     }
1319   sort_int_inv (list, 2, 1, 0, ns[1]-1);
1320
1321   for (n=0, a=0; a<ns[1]; a++)n+=list[a][0];
1322   for (best_delta=100000,best_tp=0,tp=0,fp=0,best_fp=0,a=0; a<ns[1]; a++)
1323     {
1324       tp+=list[a][0];
1325       fp+=1-list[a][0];
1326       delta=(n-(tp+fp));
1327       if (FABS(delta)<best_delta)
1328         {
1329           best_delta=delta;
1330           best_tp=tp;
1331           best_fp=fp;
1332         }
1333     }
1334   /*R[pred][real]*/
1335   tp=best_tp;
1336   fp=best_fp;
1337   fn=n-tp;
1338   tn=ns[1]-(tp+fp+fn);
1339   R[1][1]=tp;
1340   R[1][0]=fp;
1341   R[0][1]=fn;
1342   R[0][0]=tn;
1343   free_int (list, -1);
1344   return R;
1345 }
1346
1347 float** jacknife5 (Alignment*A,int **cacheIN, int *ns, int **ls, Alignment *S, char *compound, char *mode)
1348 {
1349   int seq1, ci, a,b, c, n;
1350   double score, max_score;
1351   float **R;
1352
1353
1354   int tp, tn, fn, fp;
1355   int best_tp, best_fp;
1356   int delta, best_delta;
1357   int **list;
1358   int **cache;
1359
1360   ci=name_is_in_list ( compound, S->name, S->nseq, -1);
1361   list=declare_int (A->nseq,2);
1362   R=declare_float (2, 2);
1363
1364
1365   for ( a=0; a<A->nseq; a++)
1366     {
1367       int real, res;
1368
1369       ns[0]=A->nseq-1;
1370       ns[1]=1;
1371       for (c=0,b=0; b<A->nseq; b++)
1372         if (a!=b)ls[0][c++]=b;
1373       ls[1][0]=a;
1374
1375
1376       cache=sar2cache_count1 (A, ns, ls,S, compound, mode);
1377       for (b=0; b<=26; b++)
1378         for ( c=0; c< A->len_aln; c++)
1379           cacheIN[b][c]+=cache[b][c];
1380
1381       seq1=a;
1382       real=(S->seq_al[ci][seq1]=='I')?1:0;
1383       fprintf ( stdout, ">%-10s %d ", A->name[seq1], real);
1384
1385
1386
1387       for (max_score=0,b=0; b<A->len_aln; b++)
1388         max_score+=cache[26][b];
1389
1390       for (score=0,b=0; b<A->len_aln; b++)
1391         {
1392           res=tolower (A->seq_al[seq1][b]);
1393           if ( cache[26][b]==0) continue;
1394           if ( !is_gap(res))
1395             {
1396               score+=cache[res-'a'][b];
1397             }
1398           /*fprintf ( stdout, "%c[%3d]", res,b);*/
1399         }
1400       fprintf ( stdout, " SCORE: %5d SPRED %d RATIO: %.2f \n", (int)score, a, (score*100)/max_score);
1401       list[a][0]=real;
1402
1403       if ( strstr (mode, "SIMTEST"))list[a][1]=(score*100)/max_score;
1404       else list[a][1]=(score*100)/max_score;
1405       free_int (cache, -1);
1406     }
1407
1408
1409   sort_int_inv (list, 2, 1, 0, A->nseq-1);
1410   for (n=0, a=0; a<A->nseq; a++)
1411     {
1412       n+=list[a][0];
1413     }
1414
1415   for (best_delta=100000,best_tp=0,tp=0,fp=0,best_fp=0,a=0; a<A->nseq; a++)
1416     {
1417
1418       tp+=list[a][0];
1419       fp+=1-list[a][0];
1420       delta=(n-(tp+fp));
1421       if (FABS(delta)<best_delta)
1422         {
1423           best_delta=delta;
1424           best_tp=tp;
1425           best_fp=fp;
1426         }
1427     }
1428
1429   /*R[pred][real]*/
1430   tp=best_tp;
1431   fp=best_fp;
1432   fn=n-tp;
1433   tn=A->nseq-(tp+fp+fn);
1434   R[1][1]=tp;
1435   R[1][0]=fp;
1436   R[0][1]=fn;
1437   R[0][0]=tn;
1438   free_int (list, -1);
1439
1440   return R;
1441 }
1442 float** jacknife6 (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode)
1443 {
1444   int seq1, ci, a,b, c,d,e,f, n;
1445   double score;
1446   float **R;
1447
1448
1449   int tp, tn, fn, fp;
1450   int best_tp, best_fp;
1451   int delta, best_delta;
1452   int **list;
1453
1454   ci=name_is_in_list ( compound, S->name, S->nseq, -1);
1455   list=declare_int (A->len_aln,2);
1456   R=declare_float (2, 2);
1457
1458
1459   for ( a=0; a<A->nseq; a++)
1460     {
1461       int sar, res;
1462       int **new_cache;
1463
1464       ns[0]=A->nseq-1;
1465       ns[1]=1;
1466       for (c=0,b=0; b<A->nseq; b++)
1467         if (a!=b)ls[0][c++]=b;
1468       ls[1][0]=a;
1469
1470       cache=sar2cache_proba_new (A, ns, ls,S, compound, mode);
1471
1472
1473       new_cache=declare_int (27,A->len_aln);
1474
1475       for (d=0; d< A->len_aln; d++)
1476         {
1477           int **analyze;
1478           if ( cache[26][d]==0)continue;
1479           analyze=declare_int (26, 2);
1480
1481           for ( e=0; e< ns[0]; e++)
1482             {
1483               f=ls[0][e];
1484               sar=(S->seq_al[ci][f]=='I')?1:0;
1485               res=tolower (A->seq_al[f][d]);
1486
1487               if ( res=='-') continue;
1488               analyze[res-'a'][sar]++;
1489             }
1490           for (e=0;e<26; e++)
1491             {
1492               if ( analyze[e][1]){new_cache[26][d]=1;new_cache[e][d]+=cache[e][d];}
1493               /*
1494               if ( analyze[e][0] && analyze[e][1]){new_cache[26][d]=1;new_cache[e][d]+=analyze[e][1];}
1495               else if ( analyze[e][0]){new_cache[26][d]=1;new_cache[e][d]-=analyze[e][0]*10;}
1496               else if ( analyze[e][1]){new_cache[26][d]=1;new_cache[e][d]+=analyze[e][1];}
1497               else if ( !analyze[e][0] &&!analyze[e][1]);
1498               */
1499             }
1500           free_int (analyze, -1);
1501         }
1502
1503       seq1=a;
1504       sar=(S->seq_al[ci][seq1]=='I')?1:0;
1505       fprintf ( stdout, ">%-10s %d ", A->name[seq1], sar);
1506
1507       for (score=0,b=0; b<A->len_aln; b++)
1508         {
1509           res=tolower (A->seq_al[seq1][b]);
1510           if ( cache[26][b]==0) continue;
1511           if ( !is_gap(res))
1512             {
1513               score+=new_cache[res-'a'][b];
1514             }
1515         }
1516       fprintf ( stdout, " SCORE: %5d SPRED\n", (int)score);
1517       list[seq1][0]=sar;
1518       list[seq1][1]=(int)score;
1519
1520       free_int (new_cache, -1);
1521       free_int (cache, -1);
1522     }
1523   sort_int_inv (list, 2, 1, 0, A->nseq-1);
1524   for (n=0, a=0; a<A->nseq; a++)n+=list[a][0];
1525   for (best_delta=100000,best_tp=0,tp=0,fp=0,best_fp=0,a=0; a<A->nseq; a++)
1526     {
1527       tp+=list[a][0];
1528       fp+=1-list[a][0];
1529       delta=(n-(tp+fp));
1530       if (FABS(delta)<best_delta)
1531         {
1532           best_delta=delta;
1533           best_tp=tp;
1534           best_fp=fp;
1535         }
1536     }
1537
1538   fprintf ( stderr, "\n%d %d", best_tp, best_fp);
1539   /*R[pred][real]*/
1540   tp=best_tp;
1541   fp=best_fp;
1542   fn=n-tp;
1543   tn=A->nseq-(tp+fp+fn);
1544   R[1][1]=tp;
1545   R[1][0]=fp;
1546   R[0][1]=fn;
1547   R[0][0]=tn;
1548   free_int (list, -1);
1549
1550
1551   return R;
1552 }
1553 float** cache2pred_new (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode)
1554 {
1555   int s1, seq1, ci, a,b, n;
1556   double score;
1557   float **R;
1558
1559
1560   int tp, tn, fn, fp;
1561   int best_tp, best_fp;
1562   int delta, best_delta;
1563   int **list;
1564
1565   ci=name_is_in_list ( compound, S->name, S->nseq, -1);
1566   list=declare_int ( ns[1],2);
1567   R=declare_float (2, 2);
1568
1569   for (s1=0; s1<ns[1]; s1++)
1570     {
1571       int res, real;
1572
1573       seq1=ls[1][s1];
1574       real=(S->seq_al[ci][seq1]=='I')?1:0;
1575       fprintf ( stdout, ">%-10s %d ", A->name[seq1], real);
1576       for (score=0,b=0; b<A->len_aln; b++)
1577         {
1578           res=tolower (A->seq_al[seq1][b]);
1579           if ( cache[26][b]==0) continue;
1580           if ( !is_gap(res))
1581             {
1582               score+=cache[res-'a'][b];
1583             }
1584           fprintf ( stdout, "%c", res);
1585         }
1586       fprintf ( stdout, " SCORE: %5d SPRED\n", (int)score);
1587       list[s1][0]=real;
1588       list[s1][1]=(int)score;
1589     }
1590
1591   sort_int_inv (list, 2, 1, 0, ns[1]-1);
1592
1593   for (n=0, a=0; a<ns[1]; a++)n+=list[a][0];
1594   for (best_delta=100000,best_tp=0,tp=0,fp=0,best_fp=0,a=0; a<ns[1]; a++)
1595     {
1596       tp+=list[a][0];
1597       fp+=1-list[a][0];
1598       delta=(n-(tp+fp));
1599       if (FABS(delta)<best_delta)
1600         {
1601           best_delta=delta;
1602           best_tp=tp;
1603           best_fp=fp;
1604         }
1605     }
1606
1607
1608
1609   /*R[pred][real]*/
1610   tp=best_tp;
1611   fp=best_fp;
1612   fn=n-tp;
1613   tn=ns[1]-(tp+fp+fn);
1614   R[1][1]=tp;
1615   R[1][0]=fp;
1616   R[0][1]=fn;
1617   R[0][0]=tn;
1618   free_int (list, -1);
1619
1620
1621   return R;
1622 }
1623 float** cache2pred_forbiden_res (Alignment*A,int **cache, int *ns, int **ls, Alignment *S, char *compound, char *mode)
1624 {
1625   int s1,seq1, ci, a,b, c, n;
1626   double score;
1627   float **R;
1628
1629
1630   int tp, tn, fn, fp;
1631   int best_tp, best_fp;
1632   int delta, best_delta;
1633   int **list;
1634   int **new_cache;
1635   int **mat;
1636
1637   mat=read_matrice ( "blosum62mt");
1638
1639   ci=name_is_in_list ( compound, S->name, S->nseq, -1);
1640   list=declare_int ( ns[1],2);
1641   R=declare_float (2, 2);
1642
1643   for (s1=0; s1<ns[1]; s1++)
1644     {
1645       int res, real;
1646
1647       seq1=ls[1][s1];
1648       real=(S->seq_al[ci][seq1]=='I')?1:0;
1649       fprintf ( stdout, ">%-10s %d ", A->name[seq1], real);
1650       for (score=0,b=0; b<A->len_aln; b++)
1651         {
1652           res=tolower (A->seq_al[seq1][b]);
1653           if ( cache[26][b]==0) continue;
1654           if ( !is_gap(res))
1655             {
1656               score+=cache[res-'a'][b];
1657             }
1658           fprintf ( stdout, "%c", res);
1659         }
1660       fprintf ( stdout, " SCORE: %5d SPRED\n", (int)score);
1661       list[s1][0]=real;
1662       list[s1][1]=(int)score;
1663     }
1664   new_cache=declare_int (27,A->len_aln);
1665   for (a=0; a< A->len_aln; a++)
1666     {
1667       int **analyze, real, res, d;
1668       int *res_type;
1669       int **sub;
1670       int *keep;
1671       keep=vcalloc ( 26, sizeof (int));
1672       res_type=vcalloc ( 26, sizeof (int));
1673       sub=declare_int (256, 2);
1674
1675       if ( cache[26][a]==0)continue;
1676       analyze=declare_int (26, 2);
1677       for ( b=0; b< ns[0]; b++)
1678         {
1679           seq1=ls[0][b];
1680           real=(S->seq_al[ci][seq1]=='I')?1:0;
1681           res=tolower (A->seq_al[seq1][a]);
1682
1683           if ( res=='-') continue;
1684           analyze[res-'a'][real]++;
1685         }
1686       fprintf ( stdout, "RSPRED: ");
1687       for (c=0;c<26; c++)fprintf ( stdout, "%c", c+'a');
1688       fprintf ( stdout, "\nRSPRED: ");
1689       for (c=0;c<26; c++)
1690         {
1691           if ( analyze[c][0] && analyze[c][1]){fprintf ( stdout, "1");res_type[c]='1';}
1692           else if ( analyze[c][0]){new_cache[26][a]=1;new_cache[c][a]-=analyze[c][0];fprintf ( stdout, "0");res_type[c]='0';}
1693           else if ( analyze[c][1]){new_cache[26][a]=1;new_cache[c][a]+=analyze[c][1];fprintf ( stdout, "1");res_type[c]='1';}
1694           else if ( !analyze[c][0] &&!analyze[c][1]){fprintf ( stdout, "-");res_type[c]='-';}
1695         }
1696
1697
1698       for ( c=0; c<26; c++)
1699         {
1700           for ( d=0; d<26; d++)
1701             {
1702
1703               if ( res_type[c]==res_type[d])
1704                 {
1705                   sub[res_type[c]][0]+=mat[c][d];
1706                   sub[res_type[c]][1]++;
1707                 }
1708               if ( res_type[c]!='-' && res_type[d]!='-')
1709                 {
1710                   sub['m'][0]+=mat[c][d];
1711                   sub['m'][1]++;
1712                 }
1713             }
1714         }
1715       for ( c=0; c< 256; c++)
1716         {
1717           if ( sub[c][1])fprintf ( stdout, " %c: %5.2f ", c, (float)sub[c][0]/(float)sub[c][1]);
1718         }
1719       fprintf ( stdout, " SC: %d\nRSPRED  ", cache[26][a]);
1720
1721       for ( c=0; c<26; c++)
1722         if ( res_type[c]=='1')
1723           {
1724             for (d=0; d<26; d++)
1725               if (mat[c][d]>0)keep[d]++;
1726             keep[c]=9;
1727           }
1728
1729       for (c=0; c<26; c++)
1730         {
1731           if ( keep[c]>10)fprintf ( stdout, "9");
1732           else fprintf ( stdout, "%d", keep[c]);
1733         }
1734       for ( c=0; c<26; c++)
1735         {
1736           if ( keep[c]>8)new_cache[c][a]=10;
1737           else new_cache[c][a]=-10;
1738         }
1739       fprintf ( stdout, "\n");
1740       free_int (analyze, -1);
1741       free_int (sub, -1);
1742       vfree (res_type);
1743       vfree (keep);
1744
1745     }
1746   for ( a=0; a<25; a++)
1747     for (b=a+1; b<26; b++)
1748       {
1749         int r1, r2;
1750         r1=a+'a';r2=b+'a';
1751         if ( strchr("bjoxz", r1))continue;
1752         if ( strchr("bjoxz",r2))continue;
1753
1754         if ( mat[a][b]>0 && a!=b)fprintf ( stdout, "\nMATANALYZE %c %c %d", a+'a', b+'a', mat[a][b]);
1755       }
1756
1757   for (s1=0; s1<ns[1]; s1++)
1758     {
1759       int res, real;
1760
1761       seq1=ls[1][s1];
1762       real=(S->seq_al[ci][seq1]=='I')?1:0;
1763       fprintf ( stdout, ">%-10s %d ", A->name[seq1], real);
1764       for (score=0,b=0; b<A->len_aln; b++)
1765         {
1766           res=tolower (A->seq_al[seq1][b]);
1767           if ( cache[26][b]==0) continue;
1768           if ( !is_gap(res))
1769             {
1770               score+=new_cache[res-'a'][b];
1771             }
1772           fprintf ( stdout, "%c", res);
1773         }
1774       fprintf ( stdout, " SCORE: %5d SPRED\n", (int)score);
1775       list[s1][0]=real;
1776       list[s1][1]=(int)score;
1777     }
1778   free_int (new_cache, -1);
1779   sort_int_inv (list, 2, 1, 0, ns[1]-1);
1780
1781
1782   for (n=0, a=0; a<ns[1]; a++)n+=list[a][0];
1783   for (best_delta=100000,best_tp=0,tp=0,fp=0,best_fp=0,a=0; a<ns[1]; a++)
1784     {
1785       tp+=list[a][0];
1786       fp+=1-list[a][0];
1787       delta=(n-(tp+fp));
1788       if (FABS(delta)<best_delta)
1789         {
1790           best_delta=delta;
1791           best_tp=tp;
1792           best_fp=fp;
1793         }
1794     }
1795
1796
1797
1798   /*R[pred][real]*/
1799   tp=best_tp;
1800   fp=best_fp;
1801   fn=n-tp;
1802   tn=ns[1]-(tp+fp+fn);
1803   R[1][1]=tp;
1804   R[1][0]=fp;
1805   R[0][1]=fn;
1806   R[0][0]=tn;
1807   free_int (list, -1);
1808
1809
1810   return R;
1811 }
1812
1813 int **sar2cache_proba_old ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode)
1814 {
1815   int col, s, seq,ms,mseq, res, mres, res1, n,maxn1, maxn2,maxn3, t, ci, a;
1816   float quant=0;
1817   int **list;
1818
1819   int N1msa,N1sar, N, N11, N10, N01,N00, SCORE, COL_INDEX, RES;
1820   int nfield=0;
1821   int value;
1822   float T1, T2, T3, T4;
1823   int weight_mode;
1824   int **cache;
1825   static int **sim;
1826   int sim_weight, w, sw_thr;
1827   int train_mode;
1828
1829   float zscore;
1830
1831   RES=nfield++;COL_INDEX=nfield++;N1msa=nfield++;N1sar=nfield++;N=nfield++;N11=nfield++;N10=nfield++;N01=nfield++;N00=nfield++;SCORE=nfield++;
1832   ci=name_is_in_list ( compound, S->name, S->nseq, -1);
1833   cache=declare_int (A->nseq, A->len_aln);
1834
1835   strget_param ( mode, "_FILTER1_", "0"   , "%f", &T1);
1836   strget_param ( mode, "_FILTER2_", "1000000", "%f", &T2);
1837   strget_param ( mode, "_FILTER3_", "0"   , "%f", &T3);
1838   strget_param ( mode, "_FILTER4_", "1000000", "%f", &T4);
1839   strget_param ( mode, "_SIMWEIGHT_", "1", "%d", &sim_weight);
1840   strget_param ( mode, "_SWTHR_", "30", "%d", &sw_thr);
1841   strget_param (mode, "_TRAIN_","1", "%d", &train_mode);
1842   strget_param (mode, "_ZSCORE_","0", "%f", &zscore);
1843
1844
1845
1846
1847
1848   if (sim_weight==1 && !sim) sim=aln2sim_mat(A, "idmat");
1849   for ( ms=0; ms<ns[0]; ms++)
1850     {
1851       mseq=ls[0][ms];
1852       if ( S->seq_al[ci][mseq]!='I')continue;
1853
1854       list=declare_int (A->len_aln+1, nfield);
1855       for (t=0,n=0, col=0; col< A->len_aln; col++)
1856         {
1857           int same_res;
1858
1859           mres=tolower(A->seq_al[mseq][col]);
1860           list[col][RES]=mres;
1861           list[col][COL_INDEX]=col;
1862
1863           if ( is_gap(mres))continue;
1864           for ( s=0; s<ns[0]; s++)
1865             {
1866               seq=ls[0][s];
1867               res=tolower(A->seq_al[seq][col]);
1868               if (is_gap(res))continue;
1869
1870
1871               if (sim_weight==1)
1872                 {
1873                   w=sim[seq][mseq];w=(mres==res)?100-w:w;
1874                   if (w<sw_thr)w=0;
1875                 }
1876               else
1877                 w=1;
1878
1879               if ( train_mode==4)
1880                 {
1881                   if ( S->seq_al[ci][seq]=='I')same_res=1;
1882                   else same_res=(res==mres)?1:0;
1883                 }
1884               else
1885                 same_res=(res==mres)?1:0;
1886
1887               list[col][N]+=w;
1888
1889               if (S->seq_al[ci][seq]=='I' && same_res)list[col][N11]+=w;
1890               else if (S->seq_al[ci][seq]=='I' && same_res)list[col][N10]+=w;
1891               else if (S->seq_al[ci][seq]=='O' && same_res)list[col][N01]+=w;
1892               else if (S->seq_al[ci][seq]=='O' && same_res)list[col][N00]+=w;
1893
1894               if ( S->seq_al[ci][seq]=='I')list[col][N1sar]+=w;
1895               if ( same_res)list[col][N1msa]+=w;
1896
1897             }
1898
1899           list[col][SCORE]=(int)evaluate_sar_score1 (list[col][N], list[col][N11], list[col][N1msa], list[col][N1sar]);
1900
1901         }
1902
1903       strget_param ( mode, "_MAXN1_", "5", "%d", &maxn1);
1904       strget_param ( mode, "_WEIGHT_", "1", "%d", &weight_mode);
1905       strget_param ( mode, "_QUANT_", "0.0", "%f", &quant);
1906
1907       sort_int_inv (list,nfield,SCORE,0,A->len_aln-1);
1908       if ( quant !=0)
1909         {
1910
1911           n=quantile_rank ( list,SCORE, A->len_aln,quant);
1912           sort_int (list,nfield,N1msa, 0, n-1);
1913           maxn1=MIN(n,maxn1);
1914         }
1915
1916       for (a=0; a<maxn1; a++)
1917         {
1918           col=list[a][COL_INDEX];
1919           res1=list[a][RES];
1920           value=list[a][SCORE];
1921           if ( value>T1 && value<T2){cache[mseq][col]= value;}
1922         }
1923       free_int (list, -1);
1924     }
1925
1926   /*Filter Columns*/
1927   list=declare_int (A->len_aln+1, nfield);
1928   for ( col=0; col< A->len_aln; col++)
1929     {
1930       list[col][COL_INDEX]=col;
1931       for ( s=0; s<ns[0]; s++)
1932         {
1933           seq=ls[0][s];
1934           list[col][SCORE]+=cache[seq][col];
1935         }
1936     }
1937
1938   /*Filter Columns with a score not between T2 and T3*/
1939
1940   for (col=0; col< A->len_aln; col++)
1941     if (list[col][SCORE]<T3 || list[col][SCORE]>T4)
1942       {
1943         list[col][SCORE]=0;
1944         for (s=0; s< A->nseq; s++)
1945           if (!is_gap(A->seq_al[s][col]))cache[s][col]=0;
1946       }
1947
1948   /*Keep The N Best Columns*/
1949   if ( zscore!=0)
1950     {
1951       double sum=0, sum2=0, z;
1952       int n=0;
1953       for (a=0; a< A->len_aln; a++)
1954         {
1955           if ( list[a][SCORE]>0)
1956             {
1957               sum+=list[a][SCORE];
1958               sum2+=list[a][SCORE]*list[a][SCORE];
1959               n++;
1960             }
1961         }
1962       for (a=0; a<A->len_aln; a++)
1963         {
1964           if ( list[a][SCORE]>0)
1965             {
1966               z=return_z_score (list[a][SCORE], sum, sum2,n);
1967               if ((float)z<zscore)
1968                 {
1969                   col=list[a][COL_INDEX];
1970                   for (s=0; s<A->nseq; s++)
1971                     cache [s][col]=0;
1972                 }
1973               else
1974                 {
1975                   fprintf ( stdout, "\nZSCORE: KEEP COL %d SCORE: %f SCORE: %d\n", list[a][COL_INDEX], (float)z, list[a][SCORE]);
1976                 }
1977             }
1978         }
1979     }
1980   else
1981     {
1982       sort_int_inv (list,nfield,SCORE,0,A->len_aln-1);
1983       strget_param ( mode, "_MAXN2_", "100000", "%d", &maxn2);
1984
1985       for (a=maxn2;a<A->len_aln; a++)
1986         {
1987           col=list[a][COL_INDEX];
1988           for (s=0; s<A->nseq; s++)
1989             cache [s][col]=0;
1990         }
1991     }
1992
1993   /*Get Rid of the N best Columns*/;
1994   strget_param ( mode, "_MAXN3_", "0", "%d", &maxn3);
1995
1996   for (a=0; a<maxn3;a++)
1997     {
1998       col=list[a][COL_INDEX];
1999       for (s=0; s<A->nseq; s++)
2000         cache [s][col]=0;
2001     }
2002
2003   return cache;
2004 }
2005 int aln2n_comp_col ( Alignment *A, Alignment *S, int ci)
2006 {
2007   int  res, seq,sar, col, r;
2008   int **analyze;
2009
2010   int tot=0;
2011
2012   analyze=declare_int (27, 2);
2013   for ( col=0; col< A->len_aln; col++)
2014     {
2015       int n1, n0;
2016
2017
2018       for ( n1=0, n0=0,seq=0; seq<A->nseq; seq++)
2019         {
2020           res=tolower(A->seq_al[seq][col]);
2021           sar=(S->seq_al[ci][seq]=='I')?1:0;
2022           n1+=(sar==1)?1:0;
2023           n0+=(sar==0)?1:0;
2024           if ( res=='-')continue;
2025           res-='a';
2026           analyze[res][sar]++;
2027         }
2028
2029       for (r=0; r<26; r++)
2030         {
2031           int a0,a1;
2032           a0=analyze[r][0];
2033           a1=analyze[r][1];
2034
2035
2036           if ( a1==n1 && a0<n0)
2037             {
2038               tot++;
2039             }
2040         }
2041       for ( r=0; r<26; r++)analyze[r][0]=analyze[r][1]=0;
2042     }
2043
2044   free_int (analyze, -1);
2045   return tot;
2046 }
2047 int **sar2cache_count1 ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode)
2048 {
2049   int  maxn2, res, seq,sar, ci, col,s, r;
2050   int **analyze, **list, **cache;
2051   static int **mat;
2052
2053   int a0,a1, w;
2054   if (!mat) mat=read_matrice ("blosum62mt");
2055
2056
2057   list=declare_int ( A->len_aln, 2);
2058   cache=declare_int ( 27, A->len_aln);
2059   analyze=declare_int (27, 2);
2060
2061   ci=name_is_in_list ( compound, S->name, S->nseq, -1);
2062
2063   for ( col=0; col< A->len_aln; col++)
2064     {
2065       int n1, n0;
2066
2067
2068       for ( n1=0, n0=0,s=0; s<ns[0]; s++)
2069         {
2070           seq=ls[0][s];
2071           res=tolower(A->seq_al[seq][col]);
2072           sar=(S->seq_al[ci][seq]=='I')?1:0;
2073           n1+=(sar==1)?1:0;
2074           n0+=(sar==0)?1:0;
2075           if ( res=='-')continue;
2076           res-='a';
2077
2078           analyze[res][sar]++;
2079         }
2080
2081       for (r=0; r<26; r++)
2082         {
2083
2084           a0=analyze[r][0];
2085           a1=analyze[r][1];
2086
2087           if ( strstr (mode, "SIMTEST"))
2088             {
2089               w=a1;
2090             }
2091           else if (a1 )
2092             {
2093               w=n0-a0;
2094             }
2095           else w=0;
2096
2097           cache[r][col]+=w;
2098           cache[26][col]=MAX(w, cache[26][col]);
2099         }
2100
2101       for ( r=0; r<26; r++)analyze[r][0]=analyze[r][1]=0;
2102       list[col][0]=col;
2103       list[col][1]=cache[26][col];
2104     }
2105
2106   free_int (analyze, -1);
2107
2108   sort_int_inv (list, 2, 1, 0, A->len_aln-1);
2109
2110   strget_param ( mode, "_MAXN2_", "100000", "%d", &maxn2);
2111
2112   for ( col=maxn2; col<A->len_aln; col++)
2113     for ( r=0; r<=26; r++)cache[r][list[col][0]]=0;
2114
2115   free_int (list, -1);
2116   return cache;
2117 }
2118
2119
2120 int **sar2cache_count2 ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode)
2121 {
2122   int maxn2, res, seq,sar, ci, col,s, r;
2123   int **analyze, **list, **cache, **conseq;
2124   static int **mat;
2125   int w=0;
2126   if (!mat) mat=read_matrice ("blosum62mt");
2127
2128
2129   list=declare_int ( A->len_aln, 2);
2130   cache=declare_int ( 27, A->len_aln);
2131   conseq=declare_int ( A->len_aln,3);
2132
2133   analyze=declare_int (27, 2);
2134
2135   ci=name_is_in_list ( compound, S->name, S->nseq, -1);
2136   for ( col=0; col< A->len_aln; col++)
2137     {
2138       int n1, n0;
2139
2140       for ( n1=0, n0=0,s=0; s<ns[0]; s++)
2141         {
2142           seq=ls[0][s];
2143           res=tolower(A->seq_al[seq][col]);
2144           sar=(S->seq_al[ci][seq]=='I')?1:0;
2145           n1+=(sar==1)?1:0;
2146           n0+=(sar==0)?1:0;
2147           if ( res=='-')continue;
2148           res-='a';
2149           analyze[res][sar]++;
2150         }
2151       for (r=0; r<26; r++)
2152         {
2153           int a0,a1;
2154           a0=analyze[r][0];
2155           a1=analyze[r][1];
2156           if ( a1==n1 && a0<n0)
2157             {
2158
2159               w=n0-a0;
2160               conseq[col][0]=r;
2161               conseq[col][1]=w;
2162             }
2163         }
2164       for ( r=0; r<26; r++)analyze[r][0]=analyze[r][1]=0;
2165     }
2166   free_int (analyze, -1);
2167
2168   for (s=0; s<ns[0]; s++)
2169     {
2170       int w1, w2;
2171       seq=ls[0][s];
2172       for (w1=0,w2=0,col=0; col<A->len_aln; col++)
2173         {
2174
2175           res=tolower(A->seq_al[seq][col]);
2176           if ( is_gap(res))continue;
2177           else res-='a';
2178
2179           if ( conseq[col][1] && res!=conseq[col][0])w1++;
2180           if ( conseq[col][1])w2++;
2181         }
2182       for (col=0; col<A->len_aln; col++)
2183         {
2184           res=tolower(A->seq_al[seq][col]);
2185           if ( is_gap(res))continue;
2186           else res-='a';
2187
2188           if ( conseq[col][1] && res!=conseq[col][0])conseq[col][2]+=(w2-w1);
2189         }
2190     }
2191
2192   for (col=0; col<A->len_aln; col++)
2193     {
2194       r=conseq[col][0];
2195       w=conseq[col][2];
2196
2197
2198       cache[r][col]=cache[26][col]=list[col][1]=w;
2199       list[col][0]=col;
2200     }
2201   sort_int_inv (list, 2, 1, 0, A->len_aln-1);
2202   strget_param ( mode, "_MAXN2_", "100000", "%d", &maxn2);
2203
2204   for ( col=maxn2; col<A->len_aln; col++)
2205     for ( r=0; r<=26; r++)cache[r][list[col][0]]=0;
2206
2207
2208   free_int (list, -1);
2209   return cache;
2210 }
2211
2212 int **sar2cache_count3 ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode)
2213 {
2214   int  maxn2, res, seq,sar, ci, col,s, r, a1, a0, n1, n0;
2215   int **analyze, **list, **cache;
2216   static int **mat;
2217
2218   if (!mat) mat=read_matrice ("blosum62mt");
2219
2220
2221   list=declare_int ( A->len_aln, 2);
2222   cache=declare_int ( 27, A->len_aln);
2223   analyze=declare_int (27, 2);
2224
2225   ci=name_is_in_list ( compound, S->name, S->nseq, -1);
2226
2227   for ( col=0; col< A->len_aln; col++)
2228     {
2229       double e, g;
2230       for ( n1=0, n0=0,s=0; s<ns[0]; s++)
2231         {
2232           seq=ls[0][s];
2233           res=tolower(A->seq_al[seq][col]);
2234           sar=(S->seq_al[ci][seq]=='I')?1:0;
2235           n1+=(sar==1)?1:0;
2236           n0+=(sar==0)?1:0;
2237           if ( res=='-')continue;
2238           res-='a';
2239
2240           analyze[res][sar]++;
2241         }
2242
2243       /*Gap*/
2244       for (g=0,r=0; r<A->nseq; r++)
2245         g+=is_gap(A->seq_al[r][col]);
2246       g=(100*g)/A->nseq;
2247
2248       /*enthropy
2249       for (e=0, r=0; r<26; r++)
2250         {
2251           a0=analyze[r][0];
2252           a1=analyze[r][1];
2253           t=a0+a1;
2254
2255           if (t>0)
2256             e+= t/(double)A->nseq*log(t/(double)A->nseq);
2257         }
2258       e*=-1;
2259       */
2260       e=0;
2261       if (g>10) continue;
2262       if (e>10) continue;
2263
2264       if ( strstr ( mode, "SIMTEST"))
2265         {
2266           for (r=0; r<26; r++)
2267             {
2268
2269               a0=analyze[r][0];
2270               a1=analyze[r][1];
2271
2272               if (a1)
2273                 {
2274                   cache[r][col]=a1;
2275                   cache[26][col]=MAX(cache[26][col],a1);
2276                 }
2277             }
2278         }
2279       else
2280         {
2281
2282
2283
2284           for (r=0; r<26; r++)
2285             {
2286
2287               a0=analyze[r][0];
2288               a1=analyze[r][1];
2289
2290               if (!a1 && a0)
2291                 {
2292                   cache[r][col]=a0;
2293                   cache[26][col]=MAX(cache[26][col],a0);
2294                 }
2295             }
2296         }
2297
2298       for ( r=0; r<26; r++)analyze[r][0]=analyze[r][1]=0;
2299       list[col][0]=col;
2300       list[col][1]=cache[26][col];
2301     }
2302
2303   free_int (analyze, -1);
2304
2305   sort_int_inv (list, 2, 1, 0, A->len_aln-1);
2306
2307   strget_param ( mode, "_MAXN2_", "100000", "%d", &maxn2);
2308
2309   for ( col=maxn2; col<A->len_aln; col++)
2310     for ( r=0; r<=26; r++)cache[r][list[col][0]]=0;
2311
2312   free_int (list, -1);
2313   return cache;
2314 }
2315
2316
2317 int **sar2cache_proba_new ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode)
2318 {
2319   int col, s, seq,ms,mseq, res, mres, res1, n,maxn1, maxn2,maxn3, t, ci, a,w;
2320
2321   int **list;
2322
2323   int N1msa,N1sar, N, N11, N10, N01,N00, SCORE, COL_INDEX, RES;
2324   int nfield=0;
2325   int value;
2326
2327
2328   int **cache;
2329   static int **sim;
2330   int sw_thr;
2331   float zscore;
2332
2333   RES=nfield++;COL_INDEX=nfield++;N1msa=nfield++;N1sar=nfield++;N=nfield++;N11=nfield++;N10=nfield++;N01=nfield++;N00=nfield++;SCORE=nfield++;
2334   ci=name_is_in_list ( compound, S->name, S->nseq, -1);
2335   cache=declare_int (27, A->len_aln);
2336
2337   strget_param ( mode, "_SWTHR_", "30", "%d", &sw_thr);
2338   strget_param (mode, "_ZSCORE_","0", "%f", &zscore);
2339
2340
2341   if (!sim)sim=aln2sim_mat(A, "idmat");
2342   for ( ms=0; ms<ns[0]; ms++)
2343     {
2344       mseq=ls[0][ms];
2345       if ( S->seq_al[ci][mseq]!='I')continue;
2346
2347       list=declare_int (A->len_aln+1, nfield);
2348       for (t=0,n=0, col=0; col< A->len_aln; col++)
2349         {
2350           int same_res;
2351
2352           mres=tolower(A->seq_al[mseq][col]);
2353           if ( is_gap(mres))continue;
2354
2355           list[col][RES]=mres;
2356           list[col][COL_INDEX]=col;
2357
2358           for ( s=0; s<ns[0]; s++)
2359             {
2360               seq=ls[0][s];
2361               res=tolower(A->seq_al[seq][col]);
2362               if (is_gap(res))continue;
2363               w=sim[seq][mseq];w=(mres==res)?100-w:w;
2364               if (w<sw_thr)w=0;
2365               same_res=(res==mres)?1:0;
2366
2367               list[col][N]+=w;
2368
2369               if (S->seq_al[ci][seq]=='I' && same_res)list[col][N11]+=w;
2370               else if (S->seq_al[ci][seq]=='I' && same_res)list[col][N10]+=w;
2371               else if (S->seq_al[ci][seq]=='O' && same_res)list[col][N01]+=w;
2372               else if (S->seq_al[ci][seq]=='O' && same_res)list[col][N00]+=w;
2373
2374               if ( S->seq_al[ci][seq]=='I')list[col][N1sar]+=w;
2375               if ( same_res)list[col][N1msa]+=w;
2376
2377             }
2378
2379           list[col][SCORE]=(int)evaluate_sar_score1 (list[col][N], list[col][N11], list[col][N1msa], list[col][N1sar]);
2380
2381         }
2382       strget_param ( mode, "_MAXN1_", "5", "%d", &maxn1);
2383       sort_int_inv (list,nfield,SCORE,0,A->len_aln-1);
2384       for (a=0; a<maxn1; a++)
2385         {
2386           col=list[a][COL_INDEX];
2387           res1=list[a][RES];
2388           value=list[a][SCORE];
2389
2390           if ( res1!=0)
2391             {
2392               cache[res1-'a'][col]+= value;
2393               cache[26][col]+=value;
2394             }
2395         }
2396       free_int (list, -1);
2397     }
2398
2399   /*Filter Columns*/
2400   list=declare_int (A->len_aln+1, nfield);
2401   for ( col=0; col< A->len_aln; col++)
2402     {
2403       list[col][COL_INDEX]=col;
2404       list[col][SCORE]=cache[26][col];
2405     }
2406   /*Keep The N Best Columns*/
2407   if ( zscore!=0)
2408     {
2409       double sum=0, sum2=0, z;
2410       int n=0;
2411       for (a=0; a< A->len_aln; a++)
2412         {
2413           if ( list[a][SCORE]>0)
2414             {
2415               sum+=list[a][SCORE];
2416               sum2+=list[a][SCORE]*list[a][SCORE];
2417               n++;
2418             }
2419         }
2420       for (a=0; a<A->len_aln; a++)
2421         {
2422           if ( list[a][SCORE]>0)
2423             {
2424               z=return_z_score (list[a][SCORE], sum, sum2,n);
2425               if ((float)z<zscore)
2426                 {
2427                   col=list[a][COL_INDEX];
2428                   for (s=0; s<27; s++)
2429                     cache [s][col]=0;
2430                 }
2431               else
2432                 {
2433                   fprintf ( stdout, "\nZSCORE: KEEP COL %d SCORE: %f SCORE: %d\n", list[a][COL_INDEX], (float)z, list[a][SCORE]);
2434                 }
2435             }
2436         }
2437     }
2438   else
2439     {
2440       sort_int_inv (list,nfield,SCORE,0,A->len_aln-1);
2441       strget_param ( mode, "_MAXN2_", "100000", "%d", &maxn2);
2442
2443       for (a=maxn2;a<A->len_aln; a++)
2444         {
2445           col=list[a][COL_INDEX];
2446           for (s=0; s<27; s++)
2447             cache [s][col]=0;
2448         }
2449     }
2450
2451   /*Get Rid of the N best Columns*/;
2452   strget_param ( mode, "_MAXN3_", "0", "%d", &maxn3);
2453
2454   for (a=0; a<maxn3;a++)
2455     {
2456       col=list[a][COL_INDEX];
2457       for (s=0; s<27; s++)
2458         cache [s][col]=0;
2459     }
2460   return cache;
2461 }
2462 int **sar2cache_adriana ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode)
2463 {
2464
2465   int col,maxn1, s, seq,ms,mseq, res, mres,res1, n, t, ci, a;
2466   float quant=0;
2467   int **list;
2468
2469
2470   int **cache;
2471   ci=name_is_in_list ( compound, S->name, S->nseq, -1);
2472   cache=declare_int (A->nseq, A->len_aln);
2473
2474
2475   for ( ms=0; ms<ns[0]; ms++)
2476     {
2477       mseq=ls[0][ms];
2478       if ( S->seq_al[ci][mseq]!='I')continue;
2479
2480       list=declare_int (A->len_aln+1, 5);
2481       for (t=0,n=0, col=0; col< A->len_aln; col++)
2482         {
2483           mres=tolower(A->seq_al[mseq][col]);
2484           list[col][0]=mres;
2485           list[col][1]=col;
2486
2487           if ( is_gap(mres))continue;
2488           for ( s=0; s<ns[0]; s++)
2489             {
2490               seq=ls[0][s];
2491               res=tolower(A->seq_al[seq][col]);
2492               if (is_gap(res))continue;
2493
2494               if (S->seq_al[ci][seq]=='I' && res==mres)list[col][3]++;
2495               if (res==mres)list[col][2]++;
2496             }
2497         }
2498
2499       sort_int_inv (list,5,3,0,A->len_aln-1);
2500
2501       strget_param ( mode, "_MAXN1_", "5", "%d", &maxn1);
2502       strget_param ( mode, "_QUANT_", "0.95", "%f", &quant);
2503
2504       n=quantile_rank ( list, 3, A->len_aln,quant);
2505       sort_int (list, 5, 2, 0, n-1);
2506
2507       for (a=0; a<maxn1; a++)
2508         {
2509
2510           col=list[a][1];
2511           res1=list[a][0];
2512           cache[mseq][col]=list[a][3];
2513         }
2514       free_int (list, -1);
2515
2516     }
2517   return cache;
2518 }
2519 int **sar2cache_proba2 ( Alignment *A, int *ns,int **ls, Alignment *S, char *compound, char *mode)
2520 {
2521   int col, s, seq,ms,mseq, res, mres,n,maxn1, t, ci, a,b;
2522   int COL, SCORE;
2523
2524   float quant=0;
2525   int **list;
2526
2527   float T1, T2, T3, T4;
2528
2529   int **cache;
2530   cache=declare_int ( A->nseq, A->len_aln);
2531   ci=name_is_in_list ( compound, S->name, S->nseq, -1);
2532
2533   strget_param ( mode, "_FILTER1_", "0"   , "%f", &T1);
2534   strget_param ( mode, "_FILTER2_", "1000000", "%f", &T2);
2535   strget_param ( mode, "_FILTER3_", "0"   , "%f", &T3);
2536   strget_param ( mode, "_FILTER4_", "1000000", "%f", &T4);
2537
2538   list=declare_int (A->len_aln+1,A->nseq+2);
2539   SCORE=A->nseq;
2540   COL=A->nseq+1;
2541
2542   for ( ms=0; ms<ns[0]; ms++)
2543     {
2544       mseq=ls[0][ms];
2545       if ( S->seq_al[ci][mseq]!='I')continue;
2546
2547       for (t=0,n=0, col=0; col< A->len_aln; col++)
2548         {
2549           int N11=0,N10=0,N01=0,N00=0,N1sar=0,N1msa=0,N=0;
2550
2551           mres=tolower(A->seq_al[mseq][col]);
2552           if ( is_gap(mres))continue;
2553           for ( s=0; s<ns[0]; s++)
2554             {
2555               seq=ls[0][s];
2556               res=tolower(A->seq_al[seq][col]);
2557               if (is_gap(res))continue;
2558
2559               N++;
2560               if (S->seq_al[ci][seq]=='I' && res==mres)N11++;
2561               else if (S->seq_al[ci][seq]=='I' && res!=mres)N10++;
2562               else if (S->seq_al[ci][seq]=='O' && res==mres)N01++;
2563               else if (S->seq_al[ci][seq]=='O' && res!=mres)N00++;
2564
2565               if ( S->seq_al[ci][seq]=='I')N1sar++;
2566               if ( res==mres)N1msa++;
2567             }
2568           list[col][mseq]=(int)evaluate_sar_score1 (N,N11,N1msa,N1sar);
2569           list[col][SCORE]+=list[col][mseq];
2570           list[col][COL]=col;
2571         }
2572     }
2573
2574   strget_param ( mode, "_MAXN1_", "5", "%d", &maxn1);
2575   strget_param ( mode, "_QUANT_", "0.95", "%f", &quant);
2576   sort_int_inv (list,A->nseq+2,SCORE, 0, A->len_aln-1);
2577   n=quantile_rank ( list,A->nseq, A->len_aln,quant);
2578   n=5;
2579
2580
2581   for (a=0; a<n; a++)
2582     {
2583       int value;
2584
2585       col=list[a][COL];
2586       for ( b=0; b<A->nseq; b++)
2587         {
2588           value=list[col][b];
2589           if ( value>T1 && value<T2){cache[b][col]= value;}
2590         }
2591     }
2592
2593   free_int (list, -1);
2594   return cache;
2595 }
2596
2597
2598
2599
2600 /************************************************************************************/
2601 /*                ALIGNMENT ANALYZE     : SAR                                            */
2602 /************************************************************************************/
2603 int aln2jack_group3 (Alignment *A,char *comp, int **l1, int *nl1, int **l2, int *nl2)
2604 {
2605   int **seq_list, **sar_list, nsar=0, nseq=0;
2606   int a, b, mid;
2607
2608   vsrand (0);
2609   sar_list=declare_int (A->nseq, 2);
2610   seq_list=declare_int (A->nseq, 2);
2611   for (a=0; a< A->nseq; a++)
2612     {
2613       if (comp[a]=='I')
2614         {
2615           sar_list[nsar][0]=a;
2616           sar_list[nsar][1]=rand()%100000;
2617           nsar++;
2618         }
2619       else
2620         {
2621           seq_list[nseq][0]=a;
2622           seq_list[nseq][1]=rand()%100000;
2623           nseq++;
2624         }
2625     }
2626
2627
2628   l1[0]=vcalloc (A->nseq, sizeof (int));
2629   l2[0]=vcalloc (A->nseq, sizeof (int));
2630   nl1[0]=nl2[0]=0;
2631
2632   sort_int (seq_list, 2, 1, 0,nseq-1);
2633   sort_int (sar_list, 2, 1, 0,nsar-1);
2634   mid=nsar/2;
2635   for (a=0; a<mid; a++)
2636     {
2637       l1[0][nl1[0]++]=sar_list[a][0];
2638     }
2639   for (a=0,b=mid; b<nsar; b++, a++)
2640     {
2641       l2[0][nl2[0]++]=sar_list[b][0];
2642     }
2643
2644   mid=nseq/2;
2645   for (a=0; a<mid; a++)
2646     {
2647       l1[0][nl1[0]++]=seq_list[a][0];
2648     }
2649   for (a=0,b=mid; b<nseq; b++, a++)
2650     {
2651       l2[0][nl2[0]++]=seq_list[b][0];
2652     }
2653
2654
2655   free_int (seq_list, -1);
2656   free_int (sar_list, -1);
2657   return 1;
2658 }
2659
2660 int aln2jack_group2 (Alignment *A, int seq, int **l1, int *nl1, int **l2, int *nl2)
2661 {
2662   int **list;
2663   int a, b, mid;
2664
2665
2666   list=declare_int (A->nseq, 2);
2667   l1[0]=vcalloc (A->nseq, sizeof (int));
2668   l2[0]=vcalloc (A->nseq, sizeof (int));
2669   nl1[0]=nl2[0];
2670
2671   vsrand (0);
2672   for ( a=0; a< A->nseq; a++)
2673     {
2674       list[a][0]=a;
2675       list[a][1]=rand()%100000;
2676     }
2677   sort_int (list, 2, 1, 0,A->nseq-1);
2678   mid=A->nseq/2;
2679   for (a=0; a<mid; a++)
2680     {
2681       l1[0][nl1[0]++]=list[a][0];
2682     }
2683   for (a=0,b=mid; b<A->nseq; b++, a++)
2684     {
2685       l2[0][nl2[0]++]=list[b][0];
2686     }
2687
2688   free_int (list, -1);
2689   return 1;
2690 }
2691 int aln2jack_group1 (Alignment *A, int seq, int **l1, int *nl1, int **l2, int *nl2)
2692 {
2693   int **sim;
2694   int **list;
2695   int a, mid;
2696
2697   list=declare_int ( A->nseq, 3);
2698   l1[0]=vcalloc (A->nseq, sizeof (int));
2699   l2[0]=vcalloc (A->nseq, sizeof (int));
2700   nl1[0]=nl2[0];
2701
2702   sim=aln2sim_mat (A, "idmat");
2703   for ( a=0; a< A->nseq; a++)
2704     {
2705       list[a][0]=seq;
2706       list[a][1]=a;
2707       list[a][2]=(a==seq)?100:sim[seq][a];
2708     }
2709   sort_int_inv (list, 3, 2, 0, A->nseq-1);
2710   fprintf ( stderr, "\nJacknife fromsequence %s [%d]\n", A->name[seq], seq);
2711   mid=A->nseq/2;
2712   for (a=0; a< mid; a++)
2713     l1[0][nl1[0]++]=list[a][1];
2714   for (a=mid; a<A->nseq; a++)
2715     l2[0][nl2[0]++]=list[a][1];
2716   return 1;
2717 }
2718
2719
2720 int sarset2subsarset ( Alignment *A, Alignment *S, Alignment **subA, Alignment **subS, Alignment *SUB)
2721 {
2722   Alignment *rotS, *intS;
2723   int a,b, *list, nl;
2724
2725   list=vcalloc ( SUB->nseq, sizeof (int));
2726   for (nl=0,a=0; a<SUB->nseq; a++)
2727     {
2728       b=name_is_in_list(SUB->name[a], A->name, A->nseq, 100);
2729       if ( b!=-1)list[nl++]=b;
2730     }
2731
2732   subA[0]=extract_sub_aln (A, nl, list);
2733   rotS=rotate_aln (S, NULL);
2734   intS=extract_sub_aln (rotS, nl, list);
2735
2736   subS[0]=rotate_aln (intS, NULL);
2737
2738   for ( a=0; a<S->nseq; a++) sprintf ( (subS[0])->name[a], "%s", S->name[a]);
2739
2740
2741   return 0;
2742 }
2743
2744 int ***simple_sar_analyze_vot ( Alignment *A, Alignment *SAR, char *mode)
2745 {
2746   int a, b, c, d;
2747   int res1, res2, sar1, sar2;
2748   float s;
2749   int **sim;
2750   static float ***result;
2751   static int ***iresult;
2752   if (!result)
2753     {
2754     result=declare_arrayN (3,sizeof (float),SAR->nseq, A->len_aln,3);
2755     iresult=declare_arrayN (3,sizeof (int),SAR->nseq, A->len_aln,3);
2756     }
2757
2758   sim=aln2sim_mat (A, "idmat");
2759
2760
2761   for (a=0; a<SAR->nseq; a++)
2762     for (b=0; b<A->len_aln; b++)
2763       result[a][b][0]=1;
2764
2765   for ( a=0; a< SAR->nseq; a++)
2766     for ( b=0; b<A->nseq-1; b++)
2767       for ( c=b+1; c< A->nseq; c++)
2768         for ( d=0; d<A->len_aln; d++)
2769           {
2770             res1=A->seq_al[b][d];
2771             res2=A->seq_al[c][d];
2772
2773             sar1=(SAR->seq_al[a][b]=='I')?1:0;
2774             sar2=(SAR->seq_al[a][c]=='I')?1:0;
2775
2776             s=sim[b][c];
2777
2778
2779
2780
2781             if ( sar1!=sar2 && res1!=res2)
2782               result[a][d][0]*=(1/(100-s));
2783
2784             else if ( sar1==sar2 && sar1==1 && res1==res2)
2785               result[a][d][0]*=1/s;
2786
2787
2788
2789
2790             /*
2791             else if ( sar1==sar2 && res1==res2)result[a][d][0]+=(100-s)*(100-s);
2792             else if ( sar1==sar2 && res1!=res2)result[a][d][0]-=s*s;
2793             else if ( sar1!=sar2 && res1==res2)result[a][d][0]-=(100-s)*(100-s);
2794             */
2795
2796             result[a][d][1]='a';
2797           }
2798   for ( a=0; a<SAR->nseq; a++)
2799     for ( b=0; b<A->len_aln; b++)
2800       {
2801         fprintf ( stderr, "\n%f", result[a][b][0]);
2802         iresult[a][b][0]=100*log(1-result[a][b][0]);
2803       }
2804   return iresult;
2805 }
2806 int display_simple_sar_analyze_pair_col (Alignment *A, Alignment *SAR, char *mode)
2807 {
2808   int **r;
2809   int a, b, n, do_tm;
2810
2811
2812   Alignment *rA;
2813   int *nI;
2814
2815
2816   strget_param (mode, "_TM_", "0", "%d", &do_tm);
2817   r=simple_sar_analyze_pair_col (A, SAR,mode);
2818   rA=rotate_aln (A, NULL);
2819   n=0;
2820
2821   nI=vcalloc ( SAR->nseq, sizeof (int));
2822   for (a=0; a< SAR->nseq; a++)
2823     for (b=0; b<SAR->len_aln; b++) nI[a]+=(SAR->seq_al[a][b]=='I')?1:0;
2824
2825
2826
2827   while ( r[n][0]!=-1)
2828     {
2829       if (r[n][3]>0)
2830         {fprintf ( stdout,  "COMP S: %3d %3d %s %20s %2d #\n", r[n][3],0,SAR->seq_al[r[n][0]], SAR->name[r[n][0]], nI[r[n][0]]);
2831           fprintf ( stdout, "SEQ1 S: %3d %3d %s %20s %2d #\n", r[n][3],r[n][1],rA->seq_al[r[n][1]], SAR->name[r[n][0]],nI[r[n][0]]);
2832           fprintf ( stdout, "SEQ2 S: %3d %3d %s %20s %2d #\n\n", r[n][3],r[n][2],rA->seq_al[r[n][2]], SAR->name[r[n][0]],nI[r[n][0]]);
2833         }
2834       n++;
2835     }
2836   return 0;
2837 }
2838 int display_simple_sar_analyze_col (Alignment *A, Alignment *SAR, char *mode)
2839 {
2840   int ***result, **r2, **r3, **r4, **aa;
2841   int a, b, c, n;
2842   char *cons;
2843   int threshold=20;
2844   int do_tm;
2845   strget_param (mode, "_TM_", "0", "%d", &do_tm);
2846   result=simple_sar_analyze_col (A, SAR,mode);
2847   r2=declare_int (A->len_aln*SAR->nseq, 5);
2848   r3=declare_int (A->len_aln+1, 5);
2849   r4=declare_int (A->len_aln+1, SAR->nseq+1);
2850   aa=declare_int (2, 256);
2851   cons=vcalloc (A->len_aln+1, sizeof (char));
2852   for (a=0; a<A->len_aln; a++){r3[a][0]=a;cons[a]='A';}
2853
2854
2855
2856   for (n=0,a=0; a< SAR->nseq; a++)
2857     {
2858       double sum, sum2;
2859       for (sum=0, sum2=0,b=0; b<A->len_aln; b++)
2860         {
2861           sum+=result[a][b][0];
2862           sum2+=result[a][b][0]*result[a][b][0];
2863         }
2864
2865       for (b=0; b<A->len_aln; b++, n++)
2866         {
2867           r2[n][0]=a;//compound
2868           r2[n][1]=b;//pos
2869           r2[n][2]=result[a][b][1]; //AA
2870           r2[n][3]=result[a][b][0]; //Score
2871           r2[n][4]=result[a][b][2];  //(int)10*return_z_score ((double)result[a][b][0], sum, sum2, A->len_aln); //ZScore
2872       }
2873     }
2874   sort_int (r2,5, 3, 0, n-1);//sort on Score (3rd field)
2875   for ( a=0; a< n; a++)
2876     {
2877       int comp, pos, bad;
2878
2879       comp=r2[a][0];
2880       pos=r2[a][1];
2881       fprintf ( stdout, "SEQ %5d %5d %5d %s ",r2[a][1]+1,r2[a][3], r2[a][4], (do_tm)?alnpos2hmmtop_pred (A, NULL, r2[a][1], SHORT):"x");
2882       for (c=0; c<A->nseq; c++)fprintf (stdout, "%c", A->seq_al[c][r2[a][1]]);
2883
2884
2885
2886       bad=0;
2887       for (c=0; c< A->nseq; c++)
2888         {
2889           int activity, res;
2890
2891           activity=SAR->seq_al[comp][c];
2892           res=A->seq_al[c][pos];
2893
2894           if (activity=='O')aa[0][res]++;
2895           if (activity=='I')aa[1][res]++;
2896         }
2897
2898       for (c=0; c< A->nseq; c++)
2899         {
2900           int activity, res;
2901           activity=SAR->seq_al[comp][c];
2902           res=A->seq_al[c][pos];
2903           bad+=(aa[0][res] && aa[1][res])?1:0;
2904           aa[0][res]=aa[1][res]=0;
2905         }
2906       fprintf ( stdout, " %20s %d |\nCOM %5d %5d %5d %s %s %20s %d |\n\n", SAR->name[r2[a][0]],bad,r2[a][1]+1,r2[a][3],r2[a][4], (do_tm)?alnpos2hmmtop_pred (A, NULL, r2[a][1], SHORT):"x",SAR->seq_al[r2[a][0]], SAR->name[r2[a][0]], bad);
2907
2908
2909       if (r2[a][4]>threshold)
2910         {
2911           cons[r2[a][1]]++;
2912           r3[r2[a][1]][1]++;
2913           r3[r2[a][1]][2]+=r2[a][3];
2914           r4[r2[a][1]][r2[a][0]]=1;
2915         }
2916     }
2917   sort_int (r3, 3,1,0, A->len_aln-1);
2918
2919   for (a=0; a<A->len_aln; a++)
2920     {
2921       if (r3[a][1]>0)
2922         {
2923           fprintf ( stdout, "\nPOS %4d %4d %4d %c ", r3[a][0]+1, r3[a][1], r3[a][2], cons[r3[a][0]]);
2924           for (b=0; b<SAR->nseq; b++)fprintf ( stdout, "%d", r4[r3[a][0]][b]);
2925           if (do_tm)fprintf ( stdout, " %s",alnpos2hmmtop_pred (A, NULL, r3[a][0], VERBOSE));
2926       }
2927     }
2928   for (a=0; a< A->nseq; a++)fprintf ( stdout, "\n#MSA >%s\n#MSA %s",A->name[a], A->seq_al[a]);
2929   fprintf ( stdout, "\n#MSA >cons\n#MSA %s", cons);
2930
2931   return 0;
2932 }
2933 int ***    simple_sar_predict     (Alignment *A, Alignment *SAR, char *mode)
2934 {
2935   //This function estimates the z score of every poition with every compound
2936   //The best Z-score position is then used for the prediction
2937
2938   int a, b, c, nts, pos, Rscore,Zscore;
2939   int ***r;
2940   int ***pred;
2941   int **aa;
2942
2943
2944   aa=declare_int (2,256);
2945   pred=declare_arrayN (3, sizeof (int),SAR->nseq, A->nseq, 5);
2946
2947
2948   r=simple_sar_analyze_col (A, SAR, mode);
2949   nts=SAR->len_aln; //number_trainning_sequences;
2950
2951   for (a=0; a<SAR->nseq; a++)
2952     {
2953       sort_int (r[a],4, 2, 0, A->len_aln-1);
2954
2955       pos=r[a][A->len_aln-1][3]; //Best Position
2956       Zscore=r[a][A->len_aln-1][2]; //Best Z-Score
2957       Rscore=r[a][A->len_aln-1][0]; //Best Z-Score
2958
2959
2960       for (c=0; c<nts; c++)
2961         {
2962
2963           if (SAR->seq_al[a][c]=='I')aa[1][(int)A->seq_al[c][pos]]++;//Build Positive Alphabet for Compound a
2964           if (SAR->seq_al[a][c]=='O')aa[0][(int)A->seq_al[c][pos]]++;//Build Positive Alphabet for Compound a
2965         }
2966       for (c=nts; c<A->nseq; c++)
2967         {
2968           pred[a][c][0]=pos;
2969           pred[a][c][1]=Zscore;
2970           pred[a][c][2]=Rscore;
2971           if (aa[1][(int)A->seq_al[c][pos]]>0)
2972             {
2973
2974               pred[a][c][3]=aa[1][(int)A->seq_al[c][pos]];
2975               pred[a][c][4]=aa[0][(int)A->seq_al[c][pos]];
2976             }
2977         }
2978       for (c=0; c<nts; c++)aa[0][(int)A->seq_al[c][pos]]=aa[1][(int)A->seq_al[c][pos]]=0;
2979     }
2980
2981   for ( a=nts; a< A->nseq; a++)
2982     {
2983       for ( b=0; b<SAR->nseq; b++)
2984         {
2985           fprintf ( stdout, ">%-25s %-25s Pos %3d ZScore %3d Rscore %3d Activity +: %d -: %d ", A->name [a], SAR->name[b], pred[b][a][0],pred[b][a][1], pred[b][a][2], pred[b][a][3], pred[b][a][4]);
2986           if (pred[b][a][4]==0)for (c=0; c<pred[b][a][3]; c++)fprintf ( stdout, "*");
2987           fprintf ( stdout, "\n" );
2988           for (c=0; c<A->nseq; c++)fprintf ( stdout, "%c", A->seq_al[c][ pred[b][a][0]]);
2989           fprintf ( stdout, " %s\n", SAR->name[b]);
2990           for (c=0; c<A->nseq-1; c++)fprintf ( stdout, "%c", SAR->seq_al[b][c]);
2991           fprintf ( stdout, "  %s\n", SAR->name[b]);      fprintf ( stdout, "\n");
2992
2993         }
2994     }
2995   return pred;
2996 }
2997 int *pair_seq2seq (int *iseq, char *seq1, char *seq2);
2998 int **simple_sar_analyze_pair_col ( Alignment *inA, Alignment *SAR, char *mode)
2999 {
3000
3001   int a, b, c, n, n2;
3002   int *iseq=NULL;
3003   static int **result, **fresult;
3004   int sar_mode=1;
3005   int maxgapratio=0;
3006   int nresults=10;
3007   double sum, sum2, score;
3008   Alignment *A;
3009   char aa;
3010
3011   if (!result)
3012     {
3013       result=declare_int (inA->len_aln*inA->len_aln,5);
3014
3015       fresult=declare_int (inA->len_aln*nresults*SAR->nseq, 5);
3016
3017     }
3018
3019   A=rotate_aln (inA, NULL);
3020
3021
3022   for (n2=0,a=0; a<SAR->nseq; a++)
3023     {
3024
3025       for (n=0, sum=0, sum2=0,b=0; b<A->nseq-1; b++)
3026         {
3027           for ( c=b+1; c<A->nseq; c++, n++)
3028             {
3029
3030               iseq=pair_seq2seq (iseq,A->seq_al[b], A->seq_al[c]);
3031               if ( sar_mode==1)
3032                 score=sar_vs_iseq1(SAR->seq_al[a],iseq,maxgapratio,NULL,&aa);
3033               else if (sar_mode==4)
3034                 score=sar_vs_iseq4(SAR->seq_al[a],iseq,maxgapratio,NULL,&aa);
3035               //HERE ("%d", (int)score);
3036               result[n][0]=a;//compound;
3037               result[n][1]=b; //pos1
3038               result[n][2]=c; //pos2
3039               result[n][3]=(int)score;
3040
3041               sum+=score;
3042               sum2+=score*score;
3043             }
3044         }
3045       for (b=0; b<n; b++)
3046         result[c][4]=(int)10*return_z_score ((double)result[c][3], sum, sum2,n); //ZScore
3047       sort_int (result,5,3,0,n-1);
3048       for ( b=n-nresults; b<n; b++, n2++)
3049         {
3050
3051           for ( c=0; c<5; c++)fresult[n2][c]=result[b][c];
3052         }
3053     }
3054   sort_int (fresult,5,3,0,n2-1);
3055   fresult[n2][0]=-1; //Terminator;
3056
3057   return fresult;
3058 }
3059
3060 int *pair_seq2seq (int *iseq, char *seq1, char *seq2)
3061 {
3062   static int **lu;
3063   int a;
3064   if (!lu)
3065     {
3066       int a, b, n=0;
3067       lu=declare_int (256,256);
3068       for ( n=a=0; a<256; a++)
3069         for (b=0; b<256; b++)
3070           lu[a][b]=n++;
3071     }
3072
3073   if (!iseq)iseq=vcalloc ( strlen (seq1)+1, sizeof (int));
3074   for ( a=0; a< strlen (seq1); a++)
3075     {
3076       if (is_gap(seq1[a]) || is_gap(seq2[a]))iseq[a]=-1;
3077       else iseq[a]=lu[(int)seq1[a]][(int)seq2[a]];
3078     }
3079
3080   return iseq;
3081 }
3082 int ***simple_sar_analyze_col ( Alignment *inA, Alignment *SAR, char *mode)
3083 {
3084   Alignment *A;
3085   double score=0, best_score=0;
3086   int best_pos=0;
3087   int a, b;
3088
3089   static int ***result;
3090   int **sim;
3091   char aa;
3092   int sar_mode=3;
3093   int maxgapratio=0;
3094   double sum, sum2;
3095
3096
3097   strget_param (mode, "_METHOD_", "3", "%d", &sar_mode);
3098   if (!result)
3099     result=declare_arrayN (3,sizeof (int),SAR->nseq, inA->len_aln,4);
3100
3101
3102   sim=aln2sim_mat (inA, "idmat");
3103   A=rotate_aln (inA, NULL);
3104
3105
3106   for ( a=0; a<SAR->nseq; a++)
3107     {
3108       best_pos=best_score=0;
3109       for ( sum=0, sum2=0,b=0; b<A->nseq; b++)
3110         {
3111
3112           if ( sar_mode==1)
3113             score=sar_vs_seq1(SAR->seq_al[a], A->seq_al[b],maxgapratio, sim, &aa);
3114           else if ( sar_mode==2)
3115             score=sar_vs_seq2(SAR->seq_al[a], A->seq_al[b],maxgapratio, sim, &aa);
3116           else if (sar_mode ==3)
3117             score=sar_vs_seq3(SAR->seq_al[a], A->seq_al[b],maxgapratio, sim, &aa);
3118           else if (sar_mode ==4)
3119             score=sar_vs_seq4(SAR->seq_al[a], A->seq_al[b],maxgapratio, sim, &aa);
3120           else if (sar_mode ==5)
3121             score=sar_vs_seq5(SAR->seq_al[a], A->seq_al[b],maxgapratio, sim, &aa);
3122
3123
3124           result[a][b][0]+=score*10;
3125           result[a][b][1]=aa;
3126           result[a][b][3]=b;
3127           sum+=result[a][b][0];
3128           sum2+=result[a][b][0]*result[a][b][0];
3129
3130         }
3131       for ( b=0; b< A->nseq; b++)result[a][b][2]=10*return_z_score ((double)result[a][b][0], sum, sum2, A->nseq); //Score
3132     }
3133
3134   return result;
3135
3136  }
3137 int *seq2iseq ( char *seq);
3138 double sar_vs_seq4 ( char *sar, char *seq, float gl, int **sim, char *best_aa)
3139 {
3140
3141   return sar_vs_iseq4 (sar, seq2iseq(seq), gl, sim, best_aa);
3142 }
3143 double sar_vs_seq1 ( char *sar, char *seq, float gl, int **sim, char *best_aa)
3144 {
3145
3146   return sar_vs_iseq1 (sar, seq2iseq(seq), gl, sim, best_aa);
3147 }
3148
3149 int *seq2iseq ( char *seq)
3150 {
3151  static int *iseq, clen;
3152   int a;
3153
3154   if (!iseq || clen<strlen (seq))
3155     {
3156       clen=strlen (seq);
3157       vfree (iseq);
3158       iseq=vcalloc (clen+1, sizeof (int));
3159     }
3160   for (a=0; a<strlen (seq); a++)
3161     {
3162       iseq[a]=seq[a];
3163     }
3164   return iseq;
3165 }
3166
3167 double sar_vs_iseq1 ( char *sar, int *seq, float gl, int **sim, char *best_aa)
3168 {
3169   double score=0, return_score=0;
3170   int RN,N11, Nmsa, Nsar, N, N10, N01, N00;
3171   int a, b, r, s, res, res1;
3172   double Ng=0;
3173   static int *aa, *aal, naa;
3174
3175   /*measure the E-Value for every amino acid. Returns the best one*/
3176
3177
3178
3179   N=strlen (sar);
3180   for (a=0; a<N; a++)
3181     Ng+=(is_gap(seq[a])|| seq[a]==-1);
3182   Ng/=N;
3183
3184   if (Ng>gl) return 0;
3185
3186   if (!aa)
3187     {
3188       aa=vcalloc (256*256, sizeof(int));
3189       aal=vcalloc (N, sizeof (int));
3190     }
3191   naa=0;
3192   for ( a=0; a<N; a++)
3193     {
3194       if (!aa[seq[a]])
3195         {
3196           aal[naa++]=seq[a];
3197           aa[seq[a]]=1;
3198         }
3199     }
3200
3201   best_aa[0]='-';
3202   for (a=0; a<naa; a++)
3203     {
3204
3205       res=aal[a];
3206
3207       if (res==-1 || !aa[res]);
3208       else
3209         {
3210           RN=Nmsa=Nsar=N11=N10=N01=N00=0;
3211
3212           for (b=0; b<N; b++)
3213             {
3214
3215               res1=seq[b];
3216               if (res1=='-' || res1==-1)r=0;
3217               else r=(res1==res)?1:0;
3218
3219
3220               s=(sar[b]=='I')?1:0;
3221
3222               Nmsa+=r; Nsar+=s;
3223               N11+=(r && s)?1:0;
3224               N01+=(!r &&s)?1:0;
3225               N10+=(r && !s)?1:0;
3226               N00+=(!r && !s)?1:0;
3227               RN++;
3228             }
3229
3230           if (N11)
3231             {
3232               score=evaluate_sar_score1 ( RN, N11, Nmsa, Nsar);
3233
3234             }
3235           else
3236             {
3237               score=0;
3238             }
3239
3240           if ( score>return_score)
3241             {
3242               best_aa[0]=res;
3243               return_score=score;
3244             }
3245         }
3246     }
3247   for ( a=0; a<N; a++)aa[seq[a]]=0;
3248
3249   return return_score;
3250 }
3251 double sar_vs_seq5 ( char *sar, char *seq, float gl, int **sim, char *best_aa)
3252 {
3253   double score=0;
3254   int N11, Nmsa, Nsar, N, N10, N01, N00;
3255   int a, b, r, s;
3256   double Ng=0;
3257   int *aa;
3258
3259   //measure the E-Value if all the 1AA are considered like alphabet 1*/
3260   //This function is testing ONE Compound (sar) against one column of MSA (seq) containing N sequences
3261   //Returns 0 if the active and inactive alphabet overlap
3262   N=strlen (sar);
3263   //N is the nmber of sequences
3264
3265   for (a=0; a<N; a++)
3266     Ng+=is_gap(seq[a]);
3267   Ng/=N;
3268
3269   if (Ng>gl) return 0;
3270
3271   //Identify all the AA associated with a I (Positive alphabet)
3272   aa=vcalloc ( 256, sizeof (int));
3273   for (b=0; b<N; b++)
3274     {
3275       s=(sar[b]=='I')?1:0;
3276       if (s)aa[(int)seq[b]]=1;
3277     }
3278   N11=N10=N01=N00=Nmsa=Nsar=0;
3279
3280
3281   for (b=0; b<N; b++)
3282     {
3283
3284       r=aa[(int)seq[b]];
3285       s=(sar[b]=='I')?1:0;
3286
3287       Nmsa+=r; Nsar+=s;
3288       N11+=(r && s)?1:0;
3289       N01+=(!r &&s)?1:0;
3290       N10+=(r && !s)?1:0;
3291       N00+=(!r && !s)?1:0;
3292     }
3293   if (N10>=1 || N01>=1) return 0;
3294   if (N11)
3295     {
3296       score=evaluate_sar_score1 ( N, N11, Nmsa, Nsar);
3297     }
3298   else score=0;
3299
3300   vfree (aa);
3301   return score;
3302
3303 }
3304
3305
3306
3307 double sar_vs_iseq4 ( char *sar, int *seq, float gl, int **sim, char *best_aa)
3308 {
3309   int N, Ni, No;
3310   int a, b,c, r, s;
3311   double Ng=0;
3312   static int **aa;
3313
3314   /*Correlation between AA conservation and Activity*/
3315
3316   N=strlen (sar);
3317   for (a=0; a<N; a++)
3318     Ng+=(is_gap(seq[a]) || seq[a]==-1)?1:0;
3319   Ng/=N;
3320   if (gl<1)Ng*=100;
3321
3322   if (Ng>gl) return 0;
3323
3324
3325   if (!aa)aa=declare_int(2,257*257);
3326   for (No=Ni=b=0; b<N; b++)
3327     {
3328
3329       s=(sar[b]=='I')?1:0;
3330       if (s){aa[1][seq[b]]++;Ni++;}
3331       else  {aa[0][seq[b]]++;No++;}
3332      }
3333   for (r=0,b=0; b<N; b++)
3334     {
3335       if (aa[1][(int) seq[b]]==Ni)
3336         {
3337           r=(No-aa[0][(int)seq [b]])*100/No;
3338           break;
3339         }
3340
3341     }
3342   for (c=0; c<N; c++)aa[0][seq[c]]=aa[1][seq[c]]=0;
3343
3344   return r;
3345 }
3346
3347 double sar_vs_seq3 ( char *sar, char *seq, float gl, int **sim, char *best_aa)
3348 {
3349   double score=0;
3350   int N11, Nmsa, Nsar, N, N10, N01, N00;
3351   int a, b, r, s;
3352   double Ng=0;
3353   int *aa;
3354
3355   //measure the E-Value if all the 1AA are considered like alphabet 1*/
3356   //This function is testing ONE Compound (sar) against one column of MSA (seq) containing N sequences
3357
3358   N=strlen (sar);
3359   //N is the nmber of sequences
3360
3361   for (a=0; a<N; a++)
3362     Ng+=is_gap(seq[a]);
3363   Ng/=N;
3364
3365   if (Ng>gl) return 0;
3366
3367   //Identify all the AA associated with a I (Positive alphabet)
3368   aa=vcalloc ( 256, sizeof (int));
3369   for (b=0; b<N; b++)
3370     {
3371       s=(sar[b]=='I')?1:0;
3372       if (s)aa[(int)seq[b]]=1;
3373     }
3374   N11=N10=N01=N00=Nmsa=Nsar=0;
3375
3376
3377   for (b=0; b<N; b++)
3378     {
3379
3380       r=aa[(int)seq[b]];
3381       s=(sar[b]=='I')?1:0;
3382
3383       Nmsa+=r; Nsar+=s;
3384       N11+=(r && s)?1:0;
3385       N01+=(!r &&s)?1:0;
3386       N10+=(r && !s)?1:0;
3387       N00+=(!r && !s)?1:0;
3388     }
3389
3390   if (N11)
3391     {
3392       score=evaluate_sar_score1 ( N, N11, Nmsa, Nsar);
3393     }
3394   else score=0;
3395
3396   vfree (aa);
3397   return score;
3398
3399 }
3400
3401 double sar_vs_seq2 ( char *sar, char *seq, float gl, int **sim_mat, char *best_aa)
3402 {
3403   double score=0, return_score=0;
3404   int L,N11, Nmsa, Nsar,N10, N01, N;
3405   int a, b,c,d, r1, s1,r2, s2, res;
3406   double Ng=0;
3407   int sim, diff, w;
3408   char string[5];
3409
3410   /*Weighted E-Value Similarity*/
3411   L=strlen (sar);
3412   for (a=0; a<L; a++)
3413     Ng+=is_gap(seq[a]);
3414   Ng/=L;
3415
3416   if (Ng>gl) return 0;
3417   for (a=0; a<26; a++)
3418     {
3419
3420       N=Nmsa=Nsar=N11=N10=N01=0;
3421       res='a'+a;
3422       for (d=0,b=0; b<L; b++)d+=((tolower(seq[b]))==res)?1:0;
3423       if ( d==0) continue;
3424
3425       for (b=0; b<L; b++)
3426         {
3427           r1=(tolower(seq[b])==res)?1:0;
3428           s1=(sar[b]=='I')?1:0;
3429           for ( c=0; c<L; c++)
3430             {
3431               r2=(tolower(seq[c])==res)?1:0;
3432               s2=(sar[c]=='I')?1:0;
3433
3434               sprintf ( string, "%d%d%d%d", r1,s1, r2, s2);
3435               sim= sim_mat[b][c]/10;
3436               diff=10-sim;
3437
3438               if (strm (string, "0000"))      {w=diff;N+=2*w;}
3439               else if ( strm (string, "0011")){w=sim ;N+=2*w ; N11+=w  ;N10+=0   ;N01+=w   ;Nmsa+=w   ;Nsar+=w;}
3440               else if ( strm (string, "1010")){w=diff;N+=2*w ; N11+=0  ;N10+=2*w ;N01+=0   ;Nmsa+=2*w ;Nsar+=0;}
3441               else if ( strm (string, "0101")){w=diff;N+=2*w;  N11+=0  ;N10+=0   ;N01+=2*w ;Nmsa+=0   ;Nsar+=2*w;}
3442               else if ( strm (string, "1111")){w=diff;N+=2*w;  N11+=2*w;N10+=0   ;N01+=0   ;Nmsa+=2*w ;Nsar+=2*w;}
3443               else if ( strm (string, "1001")){w=sim; N+=2*w;  N11+=0  ;N10+=w   ;N01+=w   ;Nmsa+=w;Nsar+=w;}
3444               else if ( strm (string, "0110")){w=sim; N+=2*w;  N11+=0  ;N10+=w   ;N01+=w   ;Nmsa+=w;Nsar+=w;}
3445             }
3446         }
3447       if (N11)
3448         {
3449
3450           score=evaluate_sar_score1 ( N, N11, Nmsa, Nsar);
3451         }
3452       return_score=MAX(return_score, score);
3453     }
3454   if ( return_score <0)fprintf ( stderr, "\n%.2f", return_score);
3455   return return_score;
3456 }
3457
3458 float get_sar_sim (char *seq1, char *seq2)
3459 {
3460   int a, l, s, r;
3461   int n11=0, n10=0, n01=0, n00=0;
3462
3463
3464   l=strlen (seq1);
3465   for ( a=0; a<l; a++)
3466     {
3467       s=(seq1[a]=='O')?0:1;
3468       r=(seq2[a]=='O')?0:1;
3469
3470       n00+=(!s && !r)?1:0;
3471       n11+=(s && r)?1:0;
3472       n01+=(!s && r)?1:0;
3473       n10+=(s && !r)?1:0;
3474     }
3475   if ( n11==0) return 0;
3476   else return ((float)(n11)*100)/(float)(n11+n10+n01);
3477 }
3478 float get_sar_sim2 (char *seq1, char *seq2)
3479 {
3480   int a, l, s, r;
3481   int n11=0, n10=0, n01=0, n00=0, Ns=0, Nr=0;
3482   float score1;
3483
3484   l=strlen (seq1);
3485   for ( a=0; a<l; a++)
3486     {
3487       s=(seq1[a]=='O')?0:1;
3488       r=(seq2[a]=='O')?0:1;
3489
3490       Ns+=s;
3491       Nr+=r;
3492
3493       n00+=(!s && !r)?1:0;
3494       n11+=(s && r)?1:0;
3495       n01+=(!s && r)?1:0;
3496       n10+=(s && !r)?1:0;
3497     }
3498   if ( n11==0) return 0;
3499
3500
3501
3502
3503   score1=evaluate_sar_score1 (l, n11, Ns, Nr);
3504
3505   return score1;
3506 }
3507 float sar_aln2cor (Alignment *A);
3508 int   sar_aln2ev  (Alignment *A);
3509
3510 int sarseq2ev ( char *s1, char *s2, int MODE);
3511 char*  sarseq2anti_sarseq (char *seq_in, char *seq_out);
3512 Sequence * compare_sar_sequence( Sequence *S1, Sequence *S2, int depth)
3513 {
3514   int a, b, c, d, e, f;
3515   char *s1, *s2, *s3,*s4,*s5, *n1, *n2, *n3, *n4, *n5;
3516   int max=0;;
3517   int max_depth=5;
3518   int score;
3519   char *a1, *a0, *a2, *a3, *a4;
3520   Alignment *A;
3521
3522   if ( depth>max_depth)
3523     {
3524       printf_exit (EXIT_FAILURE, stderr,"maximum depth: %d", max_depth);
3525     }
3526   if ( depth==0) depth=2;
3527   A=declare_aln2 (strlen (S1->seq[0]),depth);
3528   a0=A->seq_al[0];
3529   a1=A->seq_al[1];
3530   A->len_aln=strlen (S1->seq[0]);
3531
3532   for (a=0; a< S1->nseq; a++)
3533     for ( b=0; b<S2->nseq; b++)
3534       {
3535         A->nseq=2;
3536         sprintf (a0, "%s", S1->seq[a]);
3537         sprintf (a1, "%s", S2->seq[b]);
3538
3539         if ( strlen (a0)!=strlen (a1))
3540           {
3541             add_warning (stderr, "WARNING %s (%d) and %s (%d) do not have the same length", S1->name[a], strlen (S1->seq[a]), S2->name[b], strlen (S2->seq[b]));
3542             exit (0);
3543           }
3544
3545         fprintf ( stdout, ">2 %15s %15s CORR: %.3f EVAL: %5d\n",S1->name[a], S2->name[b],  sar_aln2cor (A), sar_aln2ev (A));
3546         sarseq2anti_sarseq (S1->seq[a],a0);
3547         fprintf ( stdout, ">2 %15s %15s ANTI: %.3f EVAL: %5d\n", S1->name[a], S2->name[b], sar_aln2cor (A), sar_aln2ev (A));
3548         if ( depth >=3)
3549           {
3550             A->nseq=3;
3551             a2=A->seq_al[2];
3552             for (c=b+1; c<S2->nseq; c++)
3553               {
3554                 sprintf (a0, "%s", S1->seq[a]);
3555                 sprintf (a1, "%s", S2->seq[b]);
3556                 sprintf (a2, "%s", S2->seq[c]);
3557                 fprintf ( stdout, ">2 %15s %15s %15s CORR: %.3f EVAL: %5d\n",S1->name[a], S2->name[b],S2->name[c],  sar_aln2cor (A), sar_aln2ev (A));
3558                 sarseq2anti_sarseq (S1->seq[a],a0);
3559                 fprintf ( stdout, ">2 %15s %15s ANTI: %.3f EVAL: %5d\n", S1->name[a], S2->name[b], S2->name[c],sar_aln2cor (A), sar_aln2ev (A));
3560               }
3561             if ( depth>=4)
3562               {
3563                 A->nseq=4;
3564                 a3=A->seq_al[2];
3565                 for (d=c+1; d<S2->nseq; d++)
3566                   {
3567                     sprintf (a0, "%s", S1->seq[a]);
3568                     sprintf (a1, "%s", S2->seq[b]);
3569                     sprintf (a2, "%s", S2->seq[c]);
3570                     sprintf (a3, "%s", S2->seq[d]);
3571
3572                     fprintf ( stdout, ">2 %15s %15s %15s %15s CORR: %.3f EVAL: %5d\n",S1->name[a], S2->name[b],S2->name[c],S2->name[d],  sar_aln2cor (A), sar_aln2ev (A));
3573                     sarseq2anti_sarseq (S1->seq[a],a0);
3574                     fprintf ( stdout, ">2 %15s %15s %15s %15s ANTI: %.3f EVAL: %5d\n", S1->name[a], S2->name[b], S2->name[c],S2->name[d],sar_aln2cor (A), sar_aln2ev (A));
3575                   }
3576                 if (depth>=5)
3577                   {
3578                     A->nseq=5;
3579                     a4=A->seq_al[3];
3580                     for (e=d+1; e<S2->nseq; e++)
3581                       {
3582                         sprintf (a0, "%s", S1->seq[a]);
3583                         sprintf (a1, "%s", S2->seq[b]);
3584                         sprintf (a2, "%s", S2->seq[c]);
3585                         sprintf (a3, "%s", S2->seq[d]);
3586                         sprintf (a4, "%s", S2->seq[d]);
3587                       }
3588                   }
3589               }
3590           }
3591       }
3592
3593   return S1;
3594 }
3595 char*  sarseq2anti_sarseq (char *seq_in, char *seq_out)
3596 {
3597   int a;
3598   if (!seq_out)seq_out=vcalloc (strlen (seq_in)+1, sizeof (char));
3599   for (a=0; a<strlen (seq_in); a++)seq_out[a]=(seq_in[a]=='I')?'O':'I';
3600   return seq_out;
3601 }
3602 float sar_aln2cor (Alignment *A)
3603 {
3604   float n1, n11, tot_cor;
3605   int a, b, c,n;
3606
3607   tot_cor=0;
3608   for (n=0,a=0; a<A->nseq-1; a++)
3609     for (b=a+1; b<A->nseq; b++)
3610       {
3611         for (n11=n1=0,c=0; c<A->len_aln; c++)
3612           {
3613             n11+=(A->seq_al[a][c]=='I' && A->seq_al[b][c]=='I');
3614             n1+= (A->seq_al[a][c]=='I' || A->seq_al[b][c]=='I');
3615           }
3616         tot_cor+=(n1==0)?0:n11/n1;
3617         n++;
3618       }
3619   tot_cor/=n;
3620   return tot_cor;
3621 }
3622 int sarseq_pair2ev  ( char *s1, char *s2,int mode);
3623 int sar_aln2ev (Alignment *A)
3624 {
3625   float n1, n11;
3626   int a, b, c, tot=0, n=0;
3627
3628   tot=0;
3629   for (a=0; a<A->nseq-1; a++)
3630     for (b=a+1; b<A->nseq; b++)
3631       {
3632         tot+=sarseq_pair2ev (A->seq_al[a], A->seq_al[b], 1);
3633         n++;
3634       }
3635   return tot;
3636 }
3637 int sarseq_pair2ev  ( char *s1, char *s2,int mode)
3638 {
3639   int l, t1, t2, t11,a, n1, n2, s;
3640   if ( (l=strlen (s1))!=strlen (s2))
3641     {
3642       return -1;
3643     }
3644   if (mode==2)
3645     {
3646       t1=l/2;
3647       t2=l/2;
3648       t11=l/2;
3649     }
3650   else
3651     {
3652       for (t1=t2=t11=0,a=0; a<l; a++)
3653         {
3654           if ( mode==1)//CORRELATED
3655             {
3656               n1=(s1[a]=='I')?1:0;
3657               n2=(s2[a]=='I')?1:0;
3658             }
3659           else if ( mode==0)//ANTICORRELATED
3660             {
3661               n1=(s1[a]=='I')?1:0;
3662               n2=(s2[a]=='O')?1:0;
3663             }
3664           t1+=n1;
3665           t2+=n2;
3666
3667           t11+=(n1&&n2)?1:0;
3668         }
3669     }
3670
3671
3672   s=(int)(100*evaluate_sar_score1 (l, t11, t1, t2));
3673   return s;
3674 }
3675
3676 double evaluate_sar_score1 ( int N, int n11, int n1msa, int n1sar)
3677 {
3678   double p, p1, p2;
3679   int n10, n01;
3680
3681   n10=n1msa-n11;
3682   n01=n1sar-n11;
3683
3684   if ( n11==0)return 0;
3685   //  if ( (n10)>n11 || n01>n11)return 0;
3686
3687   p1= M_chooses_Nlog (n1msa, N) + M_chooses_Nlog (n1sar-n11, N-n1msa) + M_chooses_Nlog (n11, n1msa);
3688   p2=(M_chooses_Nlog (n1msa, N)+  M_chooses_Nlog (n1sar, N));
3689   p=(p1-p2);
3690
3691   return -p;
3692
3693 }
3694 double evaluate_sar_score2 ( int N, int n11, int n1msa, int n1sar)
3695 {
3696
3697
3698   return n11-((n1msa-n11)+(n1sar-n11));
3699
3700   if ( n11<n1msa) return 0;
3701   else if ( n11<n1sar) return 0;
3702   else if ( n11==N)return 0;
3703   return n11;
3704 }
3705
3706
3707 int benchmark_sar( int value)
3708 {
3709   static int v[1000];
3710   static int a;
3711
3712   if (a==0)
3713     {
3714       for (a=0; a< 1000; a++)v[a]=0;
3715       v[2]=1;
3716       v[3]=2;
3717       v[6]=2;
3718       v[7]=1;
3719       v[8]=2;
3720       v[9]=1;
3721       v[10]=1;
3722       v[11]=1;
3723       v[12]=2;
3724       v[30]=2;
3725       v[31]=1;
3726       v[32]=2;
3727       v[33]=1;
3728       v[34]=2;
3729       v[35]=1;
3730       v[36]=1;
3731       v[37]=2;
3732       v[43]=2;
3733       v[44]=1;
3734       v[45]=2;
3735       v[73]=2;
3736       v[74]=1;
3737       v[75]=1;
3738       v[76]=2;
3739       v[80]=2;
3740       v[81]=1;
3741       v[82]=2;
3742       v[83]=1;
3743       v[85]=2;
3744       v[86]=1;
3745       v[87]=1;
3746       v[88]=2;
3747       v[89]=2;
3748       v[90]=1;
3749       v[91]=2;
3750       v[92]=1;
3751       v[93]=2;
3752       v[103]=2;
3753       v[104]=1;
3754       v[105]=1;
3755       v[106]=1;
3756       v[107]=2;
3757       v[130]=2;
3758       v[131]=1;
3759       v[132]=2;
3760       v[133]=1;
3761       v[134]=1;
3762       v[135]=1;
3763       v[136]=2;
3764       v[137]=1;
3765       v[138]=2;
3766       v[271]=2;
3767       v[272]=1;
3768       v[273]=2;
3769       v[281]=2;
3770       v[282]=1;
3771       v[283]=2;
3772       v[284]=1;
3773       v[285]=1;
3774       v[286]=1;
3775       v[287]=2;
3776       v[319]=2;
3777       v[320]=1;
3778       v[321]=1;
3779       v[322]=1;
3780       v[323]=1;
3781       v[324]=2;
3782       v[325]=1;
3783       v[326]=2;
3784       v[327]=1;
3785       v[328]=2;
3786       v[356]=2;
3787       v[357]=1;
3788       v[358]=1;
3789       v[359]=2;
3790       v[377]=2;
3791       v[378]=1;
3792       v[379]=2;
3793       v[386]=3;
3794       v[388]=2;
3795       v[389]=1;
3796       v[390]=1;
3797       v[391]=1;
3798       v[392]=2;
3799       v[393]=2;
3800       v[394]=2;
3801       v[395]=1;
3802       v[396]=1;
3803       v[397]=2;
3804       v[399]=2;
3805       v[400]=1;
3806       v[401]=2;
3807       v[414]=2;
3808       v[415]=1;
3809       v[416]=2;
3810       v[420]=2;
3811       v[421]=1;
3812       v[422]=1;
3813       v[423]=1;
3814       v[424]=2;
3815       v[425]=1;
3816       v[426]=2;
3817     }
3818   return v[value];
3819 }
3820
3821 Alignment *weight2sar (Alignment *A, Alignment *SAR, char *weight_file, int limit)
3822 {
3823   int a, b, c;
3824   int ***weight;
3825   char ***list;
3826   float score;
3827
3828   weight=vcalloc (SAR->nseq, sizeof (int**));
3829
3830
3831   list=file2list (weight_file, " ");
3832
3833   a=b=0;
3834   for (a=0; a< SAR->nseq; a++)
3835     {
3836       b=c=0;
3837       while (list[b])
3838         {
3839           if ( strm (list[b][1], SAR->name[a]) && atoi (list[b][3])>0)c++;
3840           b++;
3841         }
3842
3843       weight[a]=declare_int (c+1, 3);
3844       fprintf ( stderr, "\n%s %d", SAR->name[a], c);
3845       b=c=0;
3846       while (list[b])
3847         {
3848           if ( strm (list[b][1], SAR->name[a]) && atoi (list[b][3])>0)
3849             {
3850               weight[a][c][0]=atoi(list[b][2])-1;
3851               weight[a][c][1]=list[b][5][0];
3852               weight[a][c][2]=atoi (list[b][3]);
3853               c++;
3854             }
3855           b++;
3856         }
3857       weight[a][c][0]=-1;
3858     }
3859
3860   for (a=0; a<A->nseq; a++)
3861     {
3862       fprintf ( stdout, ">%s\n", A->name[a]);
3863       for ( b=0; b< SAR->nseq; b++)
3864         {
3865           score=seq2weighted_sar_score(A->seq_al[a], weight[b]);
3866           fprintf ( stdout, "%c", (score>limit)?'I':'O');
3867         }
3868       fprintf (stdout, "\n");
3869     }
3870   myexit (EXIT_SUCCESS);
3871   return A;
3872 }
3873
3874 Alignment *display_sar ( Alignment *A, Alignment *SAR, char *compound)
3875 {
3876   int a,c;
3877   char name[100];
3878
3879   c=name_is_in_list ( compound, SAR->name, SAR->nseq, 100);
3880   if ( c==-1)return A;
3881
3882   for ( a=0; a< A->nseq; a++)
3883     {
3884       sprintf (name, "%s", A->name[a]);
3885       sprintf ( A->name[a], "%c_%s_%s", SAR->seq_al[c][a], name,compound);
3886     }
3887   return A;
3888 }
3889 Alignment *aln2weighted_sar_score ( Alignment *A,Alignment *SAR, char *weight_file, char *compound)
3890 {
3891
3892   int a, b, c=0;
3893   int **weight;
3894
3895   int score;
3896   char reactivity;
3897   char ***list;
3898
3899
3900   if ( SAR)
3901     {
3902       c=name_is_in_list (compound, SAR->name, SAR->nseq, 100);
3903     }
3904
3905   list=file2list (weight_file, " ");
3906   a=b=0;
3907   while (list[a])
3908     {
3909       if (strm (list[a][1], compound))b++;
3910       a++;
3911     }
3912   weight=declare_int ( b+1, 3);
3913
3914
3915   a=b=0;
3916   while (list[a])
3917     {
3918       if ( !strm (list[a][1], compound) || strm ("TOTPOS", list[a][1]));
3919       else
3920         {
3921           weight[b][0]=atoi(list[a][2])-1;
3922           weight[b][1]=list[a][5][0];
3923           weight[b][2]=atoi(list[a][3]);
3924           b++;
3925         }
3926       a++;
3927     }
3928   weight[b][0]=-1;
3929   for ( a=0; a< A->nseq; a++)
3930     {
3931       score=seq2weighted_sar_score (A->seq_al[a], weight);
3932       reactivity=(!SAR || c==-1)?'U':SAR->seq_al[c][a];
3933
3934       sprintf (A->seq_comment[a], "Compound %-15s Reactivity %c SAR_SCORE %5d", compound,reactivity, (int) score);
3935
3936     }
3937   return A;
3938 }
3939
3940 float seq2weighted_sar_score ( char *seq, int **weight)
3941 {
3942   int a, p, r, w;
3943   float score=0;
3944
3945   a=0;
3946   while (weight[a][0]!=-1)
3947     {
3948       p=weight[a][0];
3949       r=weight[a][1];
3950       w=weight[a][2];
3951
3952       if ( is_gap(seq[p]));
3953       else if ( tolower(seq[p])==r)score+=w;
3954       a++;
3955     }
3956   return score;
3957   }
3958
3959 Alignment * sar2simpred (Alignment *A, Alignment *SAR, char *posfile, char *compound, int L1,int L2 )
3960 {
3961   int a, b, c, c1, c2;
3962   int **sim, **sim_ref, npred=0;
3963   float n11, n10, n01, n00;
3964   float sn, sp;
3965
3966   int tot_sim=0;
3967   int N11=1, N01=2, N10=3, NXX=4, SIM=5;
3968   float ***tot;
3969   int i1, i2;
3970
3971
3972   n11=n10=n01=n00=0;
3973   tot=declare_arrayN(3,sizeof (float), 10, 6, 2);
3974
3975   sim_ref=aln2sim_mat (A, "idmat");
3976   if (strm (posfile, "all"))
3977     sim=sim_ref;
3978   else
3979     {
3980       Alignment *B;
3981       B=copy_aln ( A,NULL);
3982       B=extract_aln3(B,posfile);
3983
3984       /*if (B->len_aln==0)L1=100;
3985       else
3986         L1=((B->len_aln-1)*100)/B->len_aln;
3987
3988       if (L1<=0)L1=100;
3989       */
3990       sim=aln2sim_mat (B, "idmat");
3991     }
3992
3993   for (a=0; a< A->nseq-1; a++)
3994     {
3995       for ( b=a+1; b< A->nseq; b++)
3996         {
3997           for ( c=0; c<SAR->nseq; c++)
3998             {
3999               if ( (strm (compound, SAR->name[c]) || strm ( compound, "all")))
4000                 {
4001                   /*if ( sim_ref[a][b]<30 || sim_ref[a][b]>60)continue;*/
4002                   i1=0; /*sim_ref[a][b]/10;if (i1==10)i1--;*/
4003
4004                   i2=sim[a][b];
4005
4006
4007                   c1=(SAR->seq_al[c][a]=='I')?1:0;
4008                   c2=(SAR->seq_al[c][b]=='I')?1:0;
4009
4010                   n11=(c1 && c2)?1:0;
4011                   n01=(!c1 && c2)?1:0;
4012                   n10=(c1 && !c2)?1:0;
4013                   n00=(!c1 && !c2)?1:0;
4014
4015                   tot[i1][N11][0]+=n11;
4016                   tot[i1][N01][0]+=n01;
4017                   tot[i1][N10][0]+=n10;
4018                   /*tot[i1][N00][0]+=n00;*/
4019                   tot[i1][NXX][0]++;
4020                   tot[i1][SIM][0]+=sim_ref[a][b];
4021
4022                   if ( i2>=L1)
4023                     {
4024                       tot[i1][N11][1]+=n11;
4025                       tot[i1][N01][1]+=n01;
4026                       tot[i1][N10][1]+=n10;
4027                       /*tot[i1][N00][1]+=n00;*/
4028                       tot[i1][NXX][1]++;
4029                       tot[i1][SIM][1]+=sim_ref[a][b];
4030                     }
4031                 }
4032             }
4033         }
4034     }
4035
4036   for (a=0; a<1; a++)
4037     {
4038       sp=(tot[a][N11][0])/(tot[a][N11][0]+tot[a][N10][0]);
4039       fprintf ( stdout, "\n%15s N11 %5d SP %.2f ",compound, (int)tot[a][N11][0],sp);
4040       sp=((tot[a][N11][1]+tot[a][N10][1])==0)?1:(tot[a][N11][1])/(tot[a][N11][1]+tot[a][N10][1]);
4041       sn=(tot[a][N11][0]==0)?1:(tot[a][N11][1]/tot[a][N11][0]);
4042       fprintf ( stdout, " N11 %5d SP %.2f SN %.2f SIM %.2f", (int)tot[a][N11][1], sp,sn, (tot[a][SIM][1]/tot[a][NXX][1]));
4043     }
4044
4045   myexit (0);
4046   sp=((n11+n01)==0)?1:n11/(n11+n01);
4047   sn=((n11+n01)==0)?1:n11/(n11+n10);
4048
4049   fprintf ( stdout, "\nLimit: %d NPRED %d AVGSIM %d SN %.2f   SP %.2f TP %d FP %d FN %d",L1, npred, tot_sim, sn, sp, (int)n11, (int)n01, (int)n10);
4050   myexit (EXIT_SUCCESS);
4051   return A;
4052 }
4053
4054 Alignment * sar2simpred2 (Alignment *A, Alignment *SAR, char *seqlist, char *posfile, char *compound, int L )
4055 {
4056   int a,b, c,c1, c2, p, s;
4057   float n11, n10, n01, n00, n, sn2, prediction,sp, n1, n0, t, entropy, Delta;
4058   int *rlist, *tlist, *pred, *npred, tsim, psim;
4059   int **sim, **sim_ref;
4060   int nr=0;
4061   int nrs;
4062   char *out;
4063   int delta_max;
4064   Alignment *B;
4065   int printall=1;
4066
4067   out=vcalloc (A->nseq+1, sizeof (char));
4068   rlist=vcalloc ( A->nseq, sizeof (int));
4069   tlist=vcalloc ( A->nseq, sizeof (int));
4070   pred=vcalloc(2, sizeof (int));
4071   npred=vcalloc(2, sizeof (int));
4072
4073   nrs=0;
4074   if ( strm (seqlist, "first"))
4075     {
4076       for ( a=0; a<SAR->nseq; a++)
4077         {
4078           if ( strm ( compound, SAR->name[a]))
4079             {
4080               for ( b=0; b<A->nseq; b++)
4081                 {
4082                   if ( SAR->seq_al[a][b]=='I')
4083                     {
4084                       fprintf ( stderr, "COMP: %s REF SEQ: %s\n", A->name[b], compound);
4085                       rlist[nrs]=b;
4086                       tlist[rlist[nrs]]=1;
4087                       nrs++;
4088                       break;
4089                     }
4090                 }
4091             }
4092         }
4093     }
4094   else if (strm (seqlist, "all"))
4095     {
4096       for ( a=0; a< A->nseq; a++)
4097         {
4098           rlist[nrs]=a;
4099           tlist[rlist[a]]=1;
4100           nrs++;
4101         }
4102     }
4103   else if ((a=name_is_in_list ( seqlist, A->name, A->nseq, 100))!=-1)
4104     {
4105       rlist[nrs]=a;
4106       tlist[rlist[nrs]]=1;
4107       nrs++;
4108     }
4109   else
4110     {
4111       Alignment *R;
4112       R=main_read_aln (seqlist, NULL);
4113       for (a=0; a<R->nseq; a++)
4114         {
4115           rlist[a]=name_is_in_list( R->name[a], A->name, A->nseq, 100);
4116           tlist[rlist[a]]=1;
4117         }
4118       free_aln (R);
4119     }
4120
4121   c=name_is_in_list ( compound, SAR->name, SAR->nseq, 100);
4122
4123   sim_ref=aln2sim_mat (A, "idmat");
4124   if (strm (posfile, "all"))
4125     {
4126       sim=sim_ref;
4127       B=A;
4128     }
4129   else
4130     {
4131       B=copy_aln ( A,NULL);
4132       B=extract_aln3(B,posfile);
4133       sim=aln2sim_mat (B, "idmat");
4134     }
4135
4136   n11=n10=n01=n00=n=n1=n0=0;
4137   delta_max=0;
4138   for (a=0; a<A->nseq; a++)
4139     {
4140       if ( tlist[a] && !strm (seqlist, "all"))
4141         out[a]=(SAR->seq_al[c][a]=='I')?'Z':'z';/*SAR->seq_al[c][a];*/
4142       else
4143         {
4144
4145           pred[0]=pred[1]=0;
4146           npred[0]=npred[1]=1;
4147           c1=(SAR->seq_al[c][a]=='I')?1:0;
4148           for (nr=0,tsim=0,psim=0,b=0; b<nrs; b++)
4149             {
4150               if ( SAR->seq_al[c][rlist[b]]=='o');
4151               else
4152                 {
4153                   c2=(SAR->seq_al[c][rlist[b]]=='I')?1:0;
4154                   nr+=c2;
4155                   s=sim[a][rlist[b]];
4156                   tsim+=sim_ref[a][rlist[b]];
4157                   psim+=sim[a][rlist[b]];
4158                   if (s>=L)
4159                     {
4160                       pred[c2]+=s;
4161                       npred[c2]++;
4162                     }
4163                 }
4164             }
4165
4166           if (c1==0)n0++;
4167           else n1++;
4168           t++;
4169
4170
4171           Delta=pred[1]-pred[0];
4172
4173           if (Delta<-delta_max){p=0;out[a]= (c1==0)?'O':'o';}
4174           else if (Delta>delta_max){p=1;out[a]=(c1==1)?'I':'i';}
4175           else {p=-1; out[a]=(c1==1)?'U':'u';}
4176
4177           if ( p==-1);
4178           else if (  p &&  c1)n11++;
4179           else if (  p && !c1)n10++;
4180           else if ( !p && !c1)n00++;
4181           else if ( !p &&  c1)n01++;
4182
4183           if (p!=-1)n++;
4184           if (printall)fprintf ( stdout, ">%-15s %d %c OVERALL_SIM:%d POSITION_SIM %d\n%s\n", B->name[a], c1, out[a],tsim/nrs,psim/nrs,B->seq_al[a]);
4185         }
4186     }
4187   sp=((n11+n10)==0)?1:n11/(n11+n10);
4188   sn2=((n1)==0)?1:n11/n1;
4189   prediction=(n11+n00)/(n1+n0);
4190   entropy=(float)(M_chooses_Nlog (nr, nrs)/M_chooses_Nlog(nrs/2, nrs));
4191
4192   fprintf ( stdout, ">%-15s Sp %.2f  Sn %.2f Pred %.2f E %.2f\n", compound,sp, sn2,prediction,entropy );
4193   fprintf ( stdout, "%s\n", out);
4194
4195   myexit (EXIT_SUCCESS);
4196   return A;
4197 }
4198 /************************************************************************************/
4199 /*                ALIGNMENT ANALYZE     : SAR FOR OR                                */
4200 /************************************************************************************/
4201
4202 void display_or_help();
4203 void display_or_help()
4204 {
4205   fprintf ( stdout, "\nor_sar options:");
4206   fprintf ( stdout, "\n_ORCL_: Command_line in a file");
4207
4208   fprintf ( stdout, "\n_ROTATE_  : rotate the sar matrix (if each entry is a compound rather than a sequence)");
4209   fprintf ( stdout, "\n_JNIT_    : number cycles of Jacknife");
4210   fprintf ( stdout, "\n_JNSEQ_   : Number of sequences picked up in alignment [0 to keep them all, 10 by default]");
4211   fprintf ( stdout, "\n_JSAR_    : Number of compounds picked up in the SAR matrix [0 to keep them all => default");
4212   fprintf ( stdout, "\n_JRSAR_   : Randomization of the SAR file between each JNIT iteration: S, C, R, SC, SR... (S: seq, C, column, R: residue");
4213   fprintf ( stdout, "\n_JRALN_   : Randomization of the ALN file between each JNIT iteration: S, C, R, SC, SR... (S: seq, C, column, R: residue");
4214
4215
4216   fprintf ( stdout, "\n_NPOS_    : Number of positions used to make the prediction [4 by default]");
4217   fprintf ( stdout, "\n_DEPTH_   : Depth of the motif degenerated alphabet   [Default=2]");
4218   fprintf ( stdout, "\n_POSFILE_ : Predefined list of positions in a file<P: pos score");
4219   fprintf ( stdout, "\n_RSAR_    : Randomization of the SAR file: S, C, R, SC, SR... (S: seq, C, column, R: residue");
4220   fprintf ( stdout, "\n_RALN_    : Randomization of the SAR file: S, C, R, SC, SR... (S: seq, C, column, R: residue");
4221
4222   fprintf ( stdout, "\n_SARTHR_  : E-value Threshold for the filtered displays in comploo mode (def: 3)");
4223
4224
4225
4226   fprintf ( stdout, "\n_MODE_    : jack, loo (leave one out), pos, sim, predict, self");
4227   fprintf ( stdout, "\n_COMPLIST_: Compound_list provided in a FASTA file. Each compound MUST correspond to the corresponding column in SAR");
4228   fprintf ( stdout, "\n_OUTALN_  : print the alignment corresponding to the list of positions");
4229   fprintf ( stdout, "\n_OUTTREE  : Make the tree (default nj) corrwesponding to _OUTALN_");
4230   fprintf ( stdout, "\n\n");
4231   exit (EXIT_SUCCESS);
4232 }
4233 Alignment * simple_sar_or(Alignment *inA, Alignment *inS, char *param);
4234 void sar2jack (Alignment *A, Alignment *S, int nseq, int sarlen);
4235
4236
4237 void sar2jack (Alignment *A, Alignment *S, int nseq, int sarlen)
4238 {
4239   A=aln2jacknife (A, nseq, 0);
4240   S=reorder_aln  (S, A->name, A->nseq);
4241   S=aln2jacknife (S, 0, sarlen);
4242 }
4243 Alignment *or_scan (Alignment *A,Alignment *S, char *pmode)
4244 {
4245   int l, a,ax,cx, b;
4246   char mode[100];
4247   int start, offset,w;
4248   int nl, *poslist;
4249
4250
4251   fprintf ( stdout, "\nPARAMETERS: %s\n", pmode);
4252   fprintf ( stderr, "\nPARAMETERS: %s\n", pmode);
4253
4254   strget_param (pmode, "_SW_", "15", "%d",&w);
4255   strget_param (pmode, "_SMINW_", "5", "%d",&start);
4256   strget_param (pmode, "_SSTART_", "5", "%d",&offset);
4257   strget_param (pmode, "_SMODE_", "single", "%s",&mode);
4258
4259
4260   l=intlen (A->len_aln);
4261   poslist=vcalloc ( A->len_aln, sizeof (int));
4262   nl=0;
4263
4264   for ( a=0; a< A->len_aln; a++)poslist[nl++]=a+1;
4265
4266   if ( strm (mode, "single"))
4267     {
4268       for ( a=0; a<A->len_aln-2; a++)
4269         {
4270           int c, gap, score;
4271           for (gap=0,c=0; c<A->nseq; c++)gap+=(A->seq_al[c][a]=='-');
4272           if ( !gap)
4273             {
4274               Alignment *B;
4275               B=extract_aln (A, a+1, a+2);
4276               B=or_sar (B, S, pmode, NO_PRINT);
4277               score=B->score_aln;
4278               free_aln (B);
4279             }
4280           else
4281             {
4282               score=0;
4283             }
4284           fprintf ( stdout, "P: %*d  S: %4d\n",l,a+1, score);
4285           //fprintf ( stderr, "P: %*d  S: %4d\n",l,a+1, score);
4286
4287         }
4288     }
4289   else if  strm (mode, "scan")
4290     {
4291
4292       for ( ax=0; ax<nl; ax++)
4293         {
4294           int best_score=0;
4295           int best_pos=0, best_w=0, best_start=0, best_end=0;
4296           a=poslist[ax];
4297           for (b=start; b<=w; b++)
4298             {
4299               Alignment *B;
4300               int pstart, pend;
4301               pstart=a-b;
4302               pend=a+b;
4303
4304               if (pstart<1 || pstart>A->len_aln)continue;
4305               if (pend<1 || pend>A->len_aln)continue;
4306               B=extract_aln (A, a-w, a+w);
4307
4308               B=or_sar (B, S,pmode, NO_PRINT);
4309
4310
4311               if (B->score_aln>=best_score){best_score=B->score_aln; best_pos=a;best_w=b;best_start=pstart; best_end=pend;}
4312               free_aln (B);
4313             }
4314           fprintf ( stdout, "P: %*d  I: %*d %*d Score: %5d L: %2d\n", l,best_pos, l,best_start+offset, l,best_end, best_score,(best_w*2)+1 );
4315         }
4316     }
4317   else if ( strm ( mode, "scan_comp"))
4318     {
4319       Alignment *NS=NULL;
4320       int *tresults;
4321       int s, n=0, p1;
4322       int  nbest;
4323       int *poscache;
4324       float sc, eval;
4325       int **index;
4326       int **poslist;
4327       int len=1;
4328       float best_sc, best_eval;
4329       int best_pos;
4330       char *best_word;
4331       int sth;
4332       int sev;
4333
4334       index=aln2pos_simple (A, A->nseq);
4335       strget_param (pmode, "_NPOS_", "1", "%d", &len);
4336       strget_param (pmode, "_STH_", "0", "%d", &sth);
4337       strget_param (pmode, "_SEV_", "0", "%d", &sev);
4338
4339       if (sev==0)for (a=0, sev=1; a<len; a++)sev*=A->len_aln;
4340
4341       tresults=vcalloc ( A->len_aln, sizeof (int));
4342       poscache=vcalloc ( A->len_aln, sizeof (int));
4343       poslist=generate_array_int_list (len, 0, A->len_aln-1,1, &n, NULL);
4344       for (s=0; s<S->len_aln; s++)
4345         {
4346
4347           char *word;
4348           NS=aln2block (S, s+1, s+2, NS);
4349           fprintf ( stderr, "\nProcess: %s ...\n", get_compound_name(s, mode));
4350
4351           for (best_sc=0,best_pos=0,best_word=NULL,a=0; a<n; a++)
4352             {
4353               for (b=0; b<len; b++)poscache[poslist[a][b]]=1;
4354               word=or_id_evaluate2 (A,NS, pmode,poscache, NO_PRINT, &sc);
4355               //if (word)HERE ("%s", word);
4356               eval=0;
4357               for (b=0; b<len; b++)poscache[poslist[a][b]]=0;
4358               if (word && sc>best_sc)
4359                 {
4360                   best_sc=sc;
4361                   best_eval=eval;
4362                   best_word=word;
4363                   best_pos=a;
4364                   nbest=1;
4365                 }
4366               else if ( word && sc>=best_sc)
4367                 {
4368                   nbest++;
4369                 }
4370               else
4371                 vfree (word);
4372             }
4373           nbest/=len;
4374           for (a=0; a<=A->nseq; a++)
4375             {
4376
4377
4378
4379               fprintf (stdout, "\n>%s PPPP: S: %s SC: %d EV: %d NBest: %d W: %s",  get_compound_name(s, mode),(a==A->nseq)?"cons":A->name[a],(int)best_sc, (int) best_eval, nbest,best_word);
4380               for ( b=0; b<len; b++)
4381                 {
4382                   if ( a==A->nseq)
4383                     {
4384                       p1=poslist[best_pos][b];
4385                       if (best_sc>sth && nbest<sev)tresults[p1]++;
4386                       p1++;
4387                     }
4388                   else
4389                     {
4390                       p1=index[a][poslist[best_pos][b]];
4391                       if (p1>0)p1++;
4392                     }
4393                   fprintf ( stdout, " %d ", p1);
4394
4395                 }
4396             }
4397
4398         }
4399       for ( p1=0; p1<A->len_aln;p1++)
4400         fprintf ( stdout, "\n>TOT_P: %d %d", p1+1, tresults[p1]);
4401     }
4402
4403   else if ( strm ( mode, "scan_comp_old"))
4404     {
4405       Alignment *NS=NULL, *BLOCK=NULL;
4406       int **results, *tresults;
4407       int s, n, p1, p2;
4408       int npos=1;
4409       int *poscache;
4410
4411       results=declare_int (A->len_aln*A->len_aln, 3);
4412       tresults=vcalloc ( A->len_aln, sizeof (int));
4413       poscache=vcalloc ( A->len_aln, sizeof (int));
4414       for (s=0; s<S->len_aln; s++)
4415         {
4416           int count;
4417           NS=aln2block (S, s+1, s+2, NS);
4418           fprintf ( stderr, "\nProcess: %s ...", get_compound_name(s, mode));
4419           for (n=0,p1=0; p1<A->len_aln-w; p1++, count ++)
4420             {
4421               if ( count == 50){fprintf ( stderr, "*");count=0;}
4422               for ( p2=p1; p2<A->len_aln-w; p2++, n++)
4423                 {
4424                   poscache[p1]=1;
4425                   poscache[p2]=1;
4426
4427                   BLOCK=alnpos2block (A,poscache,BLOCK);
4428
4429                   if ( aln2ngap(BLOCK)<=0)BLOCK=or_sar (BLOCK,NS, pmode, NO_PRINT);
4430                   else BLOCK->score_aln=0;
4431
4432                   //if ( BLOCK->score_aln>0)HERE ("P: %d %d %d", p1, p2,BLOCK->score_aln);
4433                   results[n][0]=p1;
4434                   results[n][1]=p2;
4435                   results[n][2]=BLOCK->score_aln;
4436                   poscache[p1]=poscache[p2]=0;
4437
4438                 }
4439             }
4440           sort_int_inv (results, 3, 2, 0, n-1);
4441           for (p1=0; p1<npos; p1++)
4442             {
4443               fprintf ( stdout, "\n>%s PPPP: %d %d SC: %d",  get_compound_name(s, mode), results[p1][0]+1,results[p1][1]+1, results[p1][2]);
4444               fprintf ( stderr, "\n>%s PPPP: %d %d SC: %d",  get_compound_name(s, mode), results[p1][0]+1,results[p1][1]+1, results[p1][2]);
4445               tresults[results[p1][0]]++;
4446               tresults[results[p1][1]]++;
4447             }
4448         }
4449       for ( p1=0; p1<A->len_aln;p1++)
4450         if ( tresults[p1])fprintf ( stdout, "\n>TOT_P: %d %d", p1+1, tresults[p1]);
4451     }
4452
4453   exit (EXIT_SUCCESS);
4454 }
4455
4456
4457 ORP * or_sar_compound(Alignment *A, Alignment *S, char *mode, int print);
4458 ORP* set_orp_name ( ORP* P, char *name);
4459 ORP* set_orp_offset ( ORP* P,int offset);
4460 Alignment * or_sar(Alignment *inA, Alignment *inS, char *mode, int print)
4461 {
4462   char rsar[4],raln[4], sarmode[100];
4463   int rotate, a, b, start, end;
4464   ORP **R, *ORS;
4465   Alignment *A, *S;
4466
4467   if ( mode && (strstr (mode, "help") || strstr (mode, "HELP")))
4468    display_or_help();
4469
4470   strget_param (mode, "_SARMODE_", "single", "%s", &sarmode);//single | all
4471   strget_param (mode, "_RSAR_", "NO", "%s", rsar);
4472   strget_param (mode, "_RALN_", "NO", "%s", raln);
4473   strget_param (mode, "_ROTATE_", "0", "%d", &rotate);
4474
4475   strget_param (mode, "_START_", "1", "%d", &start);start--;
4476   strget_param (mode, "_END_", "0", "%d", &end);
4477
4478   A=copy_aln (inA, NULL);
4479   S=copy_aln (inS, NULL);
4480
4481   if ( end==0)end=A->len_aln;
4482   if ( start!=0 || end!=A->len_aln)
4483     {
4484       int c;
4485       if ( start>=A->len_aln || end<=start || end >A->len_aln)
4486         printf_exit (EXIT_FAILURE, stderr, "ERROR: _START_%d and _END_%d are incompatible with the aln [len=%d][FATAL:%s]", start, end, A->len_aln, PROGRAM);
4487       for (b=0; b<A->nseq; b++)
4488         for (c=0,a=start; a<end; a++, c++)
4489           A->seq_al[b][c]=A->seq_al[b][a];
4490       A->len_aln=c;
4491       for (b=0; b<A->nseq; b++)A->seq_al[b][c]='\0';
4492     }
4493
4494
4495   S=aln2random_aln (S, rsar);
4496   A=aln2random_aln (A, raln);
4497
4498
4499
4500   if (rotate)
4501     {
4502       Alignment *rS;
4503
4504       if (S->len_aln!=A->nseq)
4505         printf_exit ( EXIT_FAILURE,stderr, "ERROR: Alignment and SAR matrix are incompatible [FATAL:%s]", PROGRAM);
4506
4507       rS=rotate_aln (S, NULL);
4508
4509       for (a=0; a< A->nseq; a++)sprintf (rS->name[a], "%s", A->name[a]);
4510       free_aln (S);
4511       S=rS;
4512     }
4513
4514   R=vcalloc ( S->len_aln+2, sizeof (ORP*));
4515   if (strm (sarmode, "all"))R[0]=or_sar_compound (A, S, mode, print);
4516   else if ( strm (sarmode, "single"))
4517     {
4518       for ( a=0; a<S->len_aln; a++)
4519         {
4520           Alignment *NS=NULL;
4521
4522           NS=aln2block (S, a+1, a+2, NS);
4523
4524           R[a]=or_sar_compound (A, NS, mode, NO_PRINT);
4525           set_orp_name ((R[a]), get_compound_name (a, mode));
4526           set_orp_offset(R[a], start);
4527           display_or_summary (R[a], mode, stdout, print);
4528
4529         }
4530     }
4531
4532
4533
4534   ORS=combine_n_predictions (R, A, S);
4535
4536   display_or_summary (ORS, mode, stdout, print);
4537   inA->score_aln=(int)(ORS->best*(float)1000);
4538
4539   free_orp_list (R);
4540   free_orp(ORS);
4541
4542   return inA;
4543  }
4544 ORP* set_orp_offset ( ORP* P,int offset)
4545 {
4546   if (!P) return NULL;
4547   else
4548     {
4549       P->offset=offset;
4550       return set_orp_offset(P->PR, offset);
4551     }
4552 }
4553
4554 ORP* set_orp_name ( ORP* P, char *name)
4555 {
4556   if (!P) return NULL;
4557   else
4558     {
4559       sprintf ( P->name, "%s", name);
4560       return set_orp_name(P->PR, name);
4561     }
4562 }
4563
4564 ORP * combine_n_predictions (ORP**R, Alignment *A, Alignment *S)
4565 {
4566   int a=0;
4567   ORP*N=NULL;
4568   while (R[a])
4569     {
4570       N=combine_2_predictions (R[a++], N, A, S);
4571     }
4572   sprintf ( N->name, "ALL");
4573   sprintf ( N->mode, "COMBINED");
4574
4575   return N;
4576 }
4577
4578 ORP *combine_2_predictions ( ORP*IN, ORP *TO,Alignment *A, Alignment *S)
4579 {
4580   int a;
4581
4582   if ( !TO)
4583     {
4584       TO=declare_or_prediction (IN->ncomp, IN->nseq, IN->len);
4585       TO->A=A;
4586       TO->S=S;
4587       TO->P=copy_aln(S, NULL);
4588       TO->offset=IN->offset;
4589       TO->ncomp=0;
4590     }
4591
4592   for (a=0; a< IN->len; a++)
4593     {
4594       TO->pos[a]+=IN->pos[a];
4595     }
4596
4597   TO->fp+=IN->fp;
4598   TO->fn+=IN->fn;
4599   TO->tp+=IN->tp;
4600   TO->tn+=IN->tn;
4601   rates2sensitivity (TO->tp, TO->tn, TO->fp, TO->fn, &(TO->sp), &(TO->sn), &(TO->sen2), &(TO->best));
4602
4603
4604   for (a=0; a<(TO->A)->nseq; a++)
4605     {
4606       (TO->P)->seq_al[a][TO->ncomp]=(IN->P)->seq_al[a][0];
4607       //(TO->S)->seq_al[a][TO->ncomp]=(IN->S)->seq_al[a][0];
4608     }
4609   TO->ncomp++;
4610   (TO->P)->len_aln=TO->ncomp;
4611   (TO->A)->score_aln=TO->best;
4612
4613   return TO;
4614 }
4615
4616
4617
4618 ORP * display_or_summary (ORP *CP, char *mode, FILE *fp, int print)
4619 {
4620   int a;
4621   char *pred;
4622   char *exp;
4623   char *motif;
4624   Alignment *A, *P, *S;
4625
4626
4627
4628   A=CP->A;
4629   P=CP->P;
4630   S=CP->S;
4631
4632
4633
4634   pred=vcalloc ( P->nseq*P->len_aln*2, sizeof (char));
4635   exp=vcalloc ( P->nseq*P->len_aln*2, sizeof (char));
4636   motif=vcalloc (CP->len+1, sizeof (char));
4637
4638
4639
4640   if (P && S)
4641     {
4642       for ( a=0; a<P->nseq; a++)
4643         {
4644           strcat (pred,P->seq_al[a]);
4645           strcat (exp, S->seq_al[a]);
4646         }
4647       CP->evalue=profile2evalue(pred, exp);
4648     }
4649   a=0;
4650   while ( CP->motif && CP->motif[a] && CP->motif[0][a][0])strcat (motif, CP->motif[0][a++]);
4651
4652   if ( print==PRINT)
4653     {
4654       fprintf (fp, "\n>%-10s Mode: %s Accuracy: %6.2f E-value: %6.2f Motif: ",CP->name,CP->mode, CP->best, CP->evalue);
4655
4656       if (motif[0])
4657         {
4658           fprintf (fp, " %s",motif);
4659         }
4660       for ( a=0; a<CP->len; a++)
4661         {
4662           if ( CP->pos[a]) fprintf ( fp, "\n>%-10s Mode: %s P: cons %3d SC: %4d", CP->name, CP->mode, a+1+CP->offset, CP->pos[a]);
4663         }
4664       fprintf ( fp, "\n");
4665     }
4666   vfree (pred); vfree(motif); vfree(exp);
4667   if (CP->PR)display_or_summary (CP->PR, mode, fp, print);
4668   return CP;
4669 }
4670
4671
4672 ORP * or_sar_compound(Alignment *A, Alignment *S, char *mode, int print)
4673 {
4674   char rmode[100];
4675   Alignment *P=NULL;
4676   ORP *PR=NULL;
4677
4678   strget_param (mode, "_MODE_", "predict", "%s", rmode);
4679
4680
4681
4682   if (strm (rmode, "jack"))P=or_jack (A, S, mode);
4683   else if (strm (rmode, "loo")) PR=or_loo (A, S, mode,NULL, print);
4684   else if (strm (rmode, "comploo")) P=or_comp_loo (A, S, mode,NULL, print);
4685   else if ( strm (rmode, "comppos")){or_comp_pos ( A, S, mode,print);exit(0); return NULL;}
4686   else if ( strm (rmode, "pos"))P=or_aln2pos_aln ( A, S, mode);
4687   else if ( strm (rmode, "predict"))P=or_predict ( A, S, mode);
4688   else if ( strm (rmode, "self"))PR=or_self_predict ( A, S, mode, NULL, print);
4689   else if ( strm (rmode, "sim"))P=or_sim ( A, S, mode);
4690
4691   else if ( strm (rmode, "test"))P=or_test ( A, S, mode);
4692   else
4693     {
4694        printf_exit (EXIT_FAILURE, stderr, "ERROR: %s is an unknown mode of or_sar  [FATAL:%s]\n", rmode,PROGRAM);
4695        return NULL;
4696     }
4697
4698   if (!PR)
4699     {
4700       PR=vcalloc (1, sizeof (ORP));
4701       PR->P=P;
4702     }
4703   return PR;
4704 }
4705 Alignment * or_test ( Alignment *inA, Alignment *inS, char *mode)
4706 {
4707   return inA;
4708 }
4709
4710
4711
4712 float or_id_evaluate  ( Alignment *A, Alignment *S, char *mode, int *pos, int print)
4713 {
4714   char *w;
4715   float score;
4716
4717   w=or_id_evaluate2(A,S, mode, pos,print, &score);
4718
4719   vfree (w);
4720   return score;
4721 }
4722 char* or_id_evaluate2  ( Alignment *A, Alignment *S, char *mode, int *pos, int print, float *rscore)
4723 {
4724   static char **words;
4725   static int *plist;
4726   char *bword;
4727
4728   int res, p,nl, w, c, s, exp, pred;
4729   int tp, tn, fp, fn;
4730   float sn, sp, sen2, best, best_score;
4731
4732   if (!A)
4733     {free_char (words, -1);
4734       vfree (plist);
4735       return NULL;
4736     }
4737   rscore[0]=0;
4738
4739   plist=pos2list (pos, A->len_aln, &nl);
4740   words=declare_char (A->nseq, nl+1);
4741   bword=vcalloc (nl+1, sizeof (char));
4742   for (p=0; p<nl; p++)
4743     {
4744       for (s=0; s< A->nseq; s++)
4745         {
4746           res=A->seq_al[s][plist[p]];
4747           if (res=='-'){or_id_evaluate2 (NULL, NULL, NULL, NULL, 0, NULL);vfree (bword);return 0;}
4748           words[s][p]=res;
4749         }
4750     }
4751
4752   for (best_score=0,w=0; w<A->nseq; w++)
4753     {
4754       tp=fp=fn=tn=0;
4755
4756       for (c=0; c<S->len_aln; c++)
4757         {
4758           for (s=0; s<A->nseq; s++)
4759             {
4760               exp=(S->seq_al[s][c]=='I')?1:0;
4761               pred=strm (words[w], words[s]);
4762               if      ( exp &&  pred)tp++;
4763               else if ( exp && !pred)fn++;
4764               else if (!exp && !pred)tn++;
4765               else if (!exp &&  pred)fp++;
4766             }
4767         }
4768       rates2sensitivity (tp, tn, fp,fn,&sp, &sn, &sen2, &best);
4769       if ( best>best_score)
4770         {
4771           best_score=best;
4772           sprintf (bword, "%s", words[w]);
4773         }
4774     }
4775   rscore[0]=(float)1000*best_score;
4776   or_id_evaluate2 (NULL, NULL, NULL, NULL, 0, NULL);
4777   return bword;
4778 }
4779
4780 float or_loo_evaluate2 ( Alignment *A, Alignment *S, char *mode, int *pos, int print)
4781 {
4782   int c, s, p, res, sar;
4783
4784   int *plist, nl;
4785   int tp, tn, fp, fn;
4786   float sn, sp, sen2, best;
4787   char **words, **positive, **negative;
4788
4789   tp=tn=fp=fn=0;
4790   plist=pos2list (pos, A->len_aln, &nl);
4791
4792   words=declare_char (A->nseq, nl+1);
4793
4794   for (p=0; p<nl; p++)
4795     {
4796       for (s=0; s< A->nseq; s++)
4797         {
4798           res=A->seq_al[s][plist[p]];
4799           if (res=='-'){vfree (plist);free_char (words, -1); return 0;}
4800           words[s][p]=res;
4801         }
4802     }
4803   positive=vcalloc ( A->nseq, sizeof (char*));
4804   negative=vcalloc ( A->nseq, sizeof (char*));
4805   for (c=0; c<S->len_aln; c++)
4806     {
4807       //Fill the match matrix
4808       for (p=0; p<nl; p++)
4809         {
4810           for (s=0; s< A->nseq; s++)
4811             {
4812               sar=S->seq_al[s][c];
4813               if (sar=='I')positive[s]=words[s];
4814               else if ( sar=='O')negative[s]=words[s];
4815             }
4816         }
4817
4818       //Evaluate the scores
4819       for (s=0; s< A->nseq; s++)
4820         {
4821           int pos=0, neg=0, pred;
4822           sar=S->seq_al[s][c];
4823           positive[s]=negative[s]=NULL;
4824
4825           if ( name_is_in_list (words[s], positive, A->nseq, nl+1)!=-1)
4826             pos=1;
4827           if ( name_is_in_list (words[s], negative, A->nseq, nl+1)!=-1)
4828             neg=1;
4829
4830           if (pos & !neg) pred=1;
4831           else pred=0;
4832
4833           if      ( pred  && sar=='I')tp++;
4834           else if (!pred  && sar=='I')fn++;
4835           else if (!pred  && sar=='O')tn++;
4836           else if ( pred  && sar=='O')fp++;
4837
4838           if ( sar=='I')positive[s]=words [s];
4839           else negative[s]=words[s];
4840         }
4841     }
4842
4843   vfree (negative); vfree (positive);
4844   vfree (plist); free_char (words, -1);
4845   rates2sensitivity (tp, tn, fp,fn,&sp, &sn, &sen2, &best);
4846
4847   return (float)1000*best;
4848 }
4849 float or_loo_evaluate ( Alignment *A, Alignment *S, char *mode, int *pos, int print)
4850 {
4851   int c, s, p, res, sar;
4852   int **matP,**matN;
4853   int *plist, nl;
4854   int tp, tn, fp, fn;
4855   float sn, sp, sen2, best;
4856
4857   tp=tn=fp=fn=0;
4858   plist=pos2list (pos, A->len_aln, &nl);
4859   matP=declare_int (nl, 256);
4860   matN=declare_int (nl, 256);
4861
4862   for (c=0; c<S->len_aln; c++)
4863     {
4864       //Fill the match matrix
4865       for (p=0; p<nl; p++)
4866         {
4867           for (s=0; s< A->nseq; s++)
4868             {
4869               res=A->seq_al[s][plist[p]];
4870               sar=S->seq_al[s][c];
4871               if (res=='-'){vfree (plist); free_int (matP, -1);free_int (matN, -1); return 0;}
4872               if (sar=='I')matP[p][res]++;
4873               if (sar=='O')matN[p][res]++;
4874             }
4875         }
4876
4877       //Evaluate the scores
4878       for (s=0; s< A->nseq; s++)
4879         {
4880           int scoreP, scoreN;
4881           int pred, valP, valN;
4882
4883           sar=S->seq_al[s][c];
4884           for (scoreN=0,scoreP=0,p=0; p<nl; p++)
4885             {
4886               res=A->seq_al[s][plist[p]];
4887
4888               valP=matP[p][res]-(sar=='I')?1:0;
4889               scoreP+=(valP>0)?1:0;
4890
4891               valN=matN[p][res]-(sar=='O')?1:0;
4892               scoreN+=(valN>0)?1:0;
4893             }
4894
4895           if ( scoreP==nl && scoreN<nl)pred=1;
4896           else pred=0;
4897
4898
4899           if      ( pred  && sar=='I')tp++;
4900           else if (!pred  && sar=='I')fn++;
4901           else if (!pred  && sar=='O')tn++;
4902           else if ( pred  && sar=='O')fp++;
4903         }
4904
4905       //reset the matrix
4906       for (p=0; p<nl; p++)
4907         {
4908           for (s=0; s< A->nseq; s++)
4909             {
4910               res=A->seq_al[s][plist[p]];
4911               sar=S->seq_al[s][c];
4912               if (sar=='I')matP[p][res]=0;
4913               else matN[p][res]=0;
4914             }
4915         }
4916     }
4917
4918   vfree (plist); free_int (matP, -1);free_int (matN, -1);
4919   rates2sensitivity (tp, tn, fp,fn,&sp, &sn, &sen2, &best);
4920
4921   return (float)1000*best;
4922 }
4923 int* or_comp_pos ( Alignment *inA, Alignment *inS, char *mode,int print)
4924 {
4925   Alignment *A=NULL, *S=NULL, *inS2=NULL;
4926   int a, b, c;
4927   int *main_pos, *pos=NULL;
4928
4929   set_sar (inA, inS, mode);
4930   main_pos=vcalloc ( inA->len_aln, sizeof (int));
4931
4932   inS2=copy_aln (inS, NULL);
4933   inS2->len_aln=1;
4934
4935
4936
4937   //Run every SAR, one at a time
4938   for ( c=0; c< inS->len_aln; c++)
4939     {
4940       int max, p;
4941
4942       fprintf ( stdout, ">%d\n", c);
4943       for (a=0; a< inS->nseq; a++)
4944         {
4945           inS2->seq_al[a][0]=inS->seq_al[a][c];
4946           inS2->seq_al[a][1]='\0';
4947         }
4948
4949       vfree (pos);
4950       free_aln (S);
4951       free_aln (A);
4952
4953       pos=vcalloc (inA->len_aln, sizeof (int));
4954       A=copy_aln (inA, NULL);
4955       S=copy_aln (inS2, NULL);
4956       set_sar (A,S, mode);
4957       pos=aln2predictive_positions (A, S, mode,PRINT);
4958
4959       for (max=0,b=0; b<A->len_aln; b++)
4960         {
4961           main_pos[b]+=pos[b];
4962           if (main_pos[b]>max)
4963             {
4964             max=main_pos[b];
4965             p=b+1;
4966           }
4967         }
4968
4969
4970       for (a=0; a<A->nseq; a++)
4971         {
4972           fprintf ( stdout, "\t");
4973           for ( b=0; b<A->len_aln; b++)
4974             if ( pos[b]) fprintf ( stdout, "%c", A->seq_al[a][b]);
4975           fprintf ( stdout, " %c\n", inS2->seq_al[a][0]);
4976         }
4977       fprintf ( stdout, "\n\tBest: %d %d\n", p, max);
4978
4979     }
4980
4981   if (print==PRINT)
4982     {
4983       for ( a=0; a<inA->len_aln; a++)fprintf ( stdout, "\nP2: cons %4d %4d [FINAL]", a+1, main_pos[a]);
4984     }
4985   return main_pos;
4986 }
4987 Alignment * or_comp_loo ( Alignment *inA, Alignment *inS, char *mode, int *pos,int print)
4988 {
4989   int a, b,c, n;
4990   char **keep, **remove;
4991   Alignment *A, *S, *P, *P1, *SEQ, *inS2;
4992   int **main_pos, *compound_pos, **comp_list;
4993   int pos_exists=0;
4994   char *comp_pred, *comp_exp;
4995   int sar_threshold;
4996
4997   strget_param (mode, "_SARTHRES_", "3", "%d", &sar_threshold);
4998
4999   if (pos)pos_exists=1;
5000   set_sar (inA, inS, mode);
5001   P=copy_aln (inS, NULL);
5002   keep=declare_char (inA->nseq, MAXNAMES);
5003   remove=declare_char (inA->nseq, MAXNAMES);
5004
5005   main_pos=declare_int ( inA->len_aln,4);
5006   comp_list=declare_int (inA->len_aln, sizeof (int*));
5007   inS2=copy_aln (inS, NULL);
5008   inS2->len_aln=1;
5009
5010
5011   comp_pred=vcalloc ( inA->nseq+1, sizeof (int));
5012   comp_exp=vcalloc ( inA->nseq+1, sizeof (int));
5013   compound_pos=NULL;
5014   //Run every SAR, one at a time
5015   for ( c=0; c< inS->len_aln; c++)
5016     {
5017       for (a=0; a< inS->nseq; a++)
5018         {
5019           inS2->seq_al[a][0]=inS->seq_al[a][c];
5020           inS2->seq_al[a][1]='\0';
5021         }
5022       vfree (compound_pos);
5023       compound_pos=vcalloc (inA->len_aln, sizeof (int));
5024       for (a=0; a<inA->nseq; a++)
5025         {
5026           char ***motifs;
5027
5028           A=copy_aln (inA, NULL);
5029           S=copy_aln (inS2, NULL);
5030
5031           for (n=0,b=0; b<A->nseq; b++)
5032             {
5033               if (b!=a)sprintf (keep[n++ ], "%s", A->name [b]);
5034             }
5035           sprintf ( remove[0], "%s", A->name[a]);
5036           reorder_aln (A,keep, A->nseq-1);
5037
5038           set_sar (A,S, mode);
5039           if (!pos_exists)
5040             {
5041               pos=aln2predictive_positions (A, S, mode,NO_PRINT);
5042             }
5043           for (b=0; b<A->len_aln; b++)
5044             {
5045               compound_pos[b]+=pos[b];
5046             }
5047
5048           motifs=compounds2motifs (A, S, pos,0, mode, NO_PRINT);
5049
5050           SEQ=copy_aln (inA, NULL);
5051           SEQ=reorder_aln (SEQ, remove, 1);
5052
5053           P1=aln2prediction (SEQ, motifs, pos);
5054           comp_pred[a]=P1->seq_al[0][0];
5055           comp_exp[a]=inS2->seq_al[a][0];
5056           P->seq_al[a][c]=P1->seq_al[0][0];
5057
5058
5059           free_aln (SEQ);
5060           free_aln (S);
5061           free_aln (A);
5062           free_aln (P1);
5063           free_arrayN( (void *)motifs, 3);
5064           if (!pos_exists)vfree (pos);
5065         }
5066       if (print==PRINT)
5067         fprintf ( stdout, ">%-15s SC: %.2f E; %.2f\n%s\n%s\n", get_compound_name(c, mode),profile2sensitivity (comp_pred, comp_exp, NULL, NULL, NULL, NULL),profile2evalue(comp_pred, comp_exp),comp_pred, comp_exp);
5068       for (b=0; b<A->len_aln; b++)
5069             {
5070               main_pos[b][2]+=compound_pos[b];
5071               if (compound_pos[b])main_pos[b][3]++;
5072               if (profile2evalue(comp_pred, comp_exp)>sar_threshold)
5073                 {
5074                   main_pos[b][0]+=compound_pos[b];
5075                   if (compound_pos[b])
5076                     {
5077                       main_pos[b][1]++;
5078                       comp_list[b][0]++;
5079                       comp_list[b]=vrealloc (comp_list[b], sizeof (int)*(comp_list[b][0]+1));
5080                       comp_list[b][comp_list[b][0]]=c;
5081                     }
5082                 }
5083             }
5084
5085     }
5086
5087   P->score_aln=(int)((float)1000*evaluate_prediction (P, inS, mode,print));
5088
5089   if (print==PRINT)
5090     {
5091       for ( a=0; a<inA->len_aln; a++)fprintf ( stdout, "\nP: cons %4d RS: %4d RC: %5d FC: %4d %4d[FINAL]", a+1, main_pos[a][2], main_pos[a][3], main_pos[a][0], main_pos[a][1]);
5092
5093       for ( a=0; a<inA->len_aln; a++)
5094         {
5095           fprintf ( stdout, "\nP: cons %4d RS: %4d RC: %5d FC: %4d %4d CLIST: ", a+1, main_pos[a][2], main_pos[a][3], main_pos[a][0], main_pos[a][1]);
5096           for ( c=1; c<=comp_list[a][0]; c++)
5097             {
5098               fprintf ( stdout, "%s ",  get_compound_name(comp_list[a][c], mode));
5099             }
5100           fprintf ( stdout, " [COMP_LIST]");
5101         }
5102     }
5103   free_int (main_pos, -1);
5104   free_int (comp_list, -1);
5105   return P;
5106 }
5107
5108 ORP* or_loo ( Alignment *inA, Alignment *inS, char *mode, int *pos,int print)
5109 {
5110   int a, b,  n;
5111   char **keep, **remove;
5112   Alignment *A, *S, *P, *P1, *SEQ;
5113
5114   int pos_exists=0;
5115   ORP *PR;
5116
5117
5118
5119   if (pos)pos_exists=1;
5120   set_sar (inA, inS, mode);
5121   PR=declare_or_prediction (inS->nseq, inA->nseq, inA->len_aln);
5122   sprintf (PR->mode, "loo ");
5123   P=copy_aln (inS, NULL);
5124
5125
5126   keep=declare_char (inA->nseq, MAXNAMES);
5127   remove=declare_char (inA->nseq, MAXNAMES);
5128
5129
5130   PR->A=inA;
5131   PR->P=P;
5132   PR->S=inS;
5133
5134   for (a=0; a<inA->nseq; a++)
5135     {
5136       char ***motifs;
5137
5138       A=copy_aln (inA, NULL);
5139       S=copy_aln (inS, NULL);
5140
5141       for (n=0,b=0; b<A->nseq; b++)
5142         {
5143           if (b!=a)sprintf (keep[n++ ], "%s", A->name [b]);
5144         }
5145       sprintf ( remove[0], "%s", A->name[a]);
5146       reorder_aln (A,keep, A->nseq-1);
5147
5148       set_sar (A,S, mode);
5149
5150       if (!pos_exists)
5151         {
5152           pos=aln2predictive_positions (A, S, mode,print);
5153
5154         }
5155
5156       for (b=0; b<A->len_aln; b++)
5157         {
5158           PR->pos[b]+=pos[b];
5159         }
5160
5161       motifs=compounds2motifs (A, S, pos,0, mode, print);
5162
5163       SEQ=copy_aln (inA, NULL);
5164       SEQ=reorder_aln (SEQ, remove, 1);
5165       SEQ->nseq=1;
5166
5167       P1=aln2prediction (SEQ, motifs, pos);
5168
5169
5170       if (print==PRINT)
5171         {
5172           fprintf ( stdout, "\n%s\nPred: %s\nReal: %s\n", P1->name[0], P1->seq_al[0], inS->seq_al[a]);
5173         }
5174       sprintf ( P->seq_al[a], "%s", P1->seq_al[0]);
5175       free_aln (P1);
5176
5177
5178       free_aln (SEQ);
5179       free_aln (S);
5180       free_aln (A);
5181
5182       free_arrayN( (void *)motifs, 3);
5183       if (!pos_exists)vfree (pos);
5184
5185     }
5186   free_char (keep, -1);
5187   free_char (remove, -1);
5188
5189
5190   PR=new_evaluate_prediction (PR, mode,print);
5191
5192
5193   PR->PR=or_self_predict(inA, inS, mode,NULL, print);
5194
5195   if (print==PRINT)for ( a=0; a<inA->len_aln; a++)fprintf ( stdout, "\nP: cons %d %d [FINAL]", a+1, PR->pos[a]);
5196
5197
5198
5199   return PR;
5200 }
5201
5202
5203
5204 Alignment * or_jack(Alignment *inA, Alignment *inS, char *mode)
5205 {
5206   int a,b;
5207   int n_cycles=100;
5208   int subnseq=10;
5209   int subsar=0;
5210   Alignment *A, *S;
5211   int *main_pos,*pos;
5212   char jrsar[10], jraln[10];
5213
5214   strget_param (mode, "_JNIT_", "100", "%d", &n_cycles);
5215   strget_param (mode, "_JNSEQ_", "10", "%d", &subnseq);
5216   strget_param (mode, "_JNAR_", "0", "%d", &subsar);
5217
5218   strget_param (mode, "_JRSAR_", "NO", "%s", jrsar);
5219   strget_param (mode, "_JRALN_", "NO", "%s", jraln);
5220
5221
5222
5223   main_pos=vcalloc ( inA->len_aln, sizeof (int));
5224   for (a=0; a<n_cycles; a++)
5225     {
5226
5227       A=copy_aln (inA, NULL);
5228       S=copy_aln (inS, NULL);
5229
5230       S=aln2random_aln (S, jrsar);
5231       A=aln2random_aln (A, jraln);
5232
5233       set_sar (A,S, mode);
5234       sar2jack (A, S,subnseq,subsar);
5235
5236
5237
5238       pos=aln2predictive_positions (A, S,mode, PRINT);
5239
5240       for (b=0; b<A->len_aln; b++)main_pos[b]+=pos[b];
5241       vfree (pos);
5242
5243     }
5244   display_pos (A, S, main_pos, mode);
5245
5246
5247   return inA;
5248 }
5249
5250 Alignment * display_pos (Alignment *A, Alignment *S, int *pos,char *mode)
5251 {
5252   Alignment *B;
5253   int a, b;
5254   int **index;
5255
5256   int intl;
5257
5258   intl=intlen (A->len_aln);
5259   index=aln2pos_simple (A, A->nseq);
5260   B=copy_aln (A,NULL);
5261   B->len_aln=0;
5262   for ( a=0; a<A->len_aln; a++)
5263     fprintf ( stdout, "\nP: cons %*d %*d S: %4d [DISPLAY_FULL_POS]", intl,a+1,intl, a+2, pos[a]);
5264   fprintf ( stdout, "\n\n");
5265   for (a=0; a<A->len_aln; a++)
5266     {
5267       if (pos[a])
5268         {
5269           for ( b=0; b<A->nseq; b++)
5270             {
5271               B->seq_al[b][B->len_aln]=A->seq_al[b][a];
5272               if (index[b][a]>0)fprintf ( stdout, "\nP: %s %d %d S: %d [DISPLAY_POS]",A->name[b], index[b][a], index[b][a]+1, pos[a]);
5273             }
5274           B->len_aln++;
5275           fprintf ( stdout, "\nP: cons %d %d S: %d [DISPLAY_POS]", a+1, a+2, pos[a]);
5276         }
5277     }
5278   fprintf ( stdout, "\n");
5279   for (a=0; a<B->nseq; a++)B->seq_al[a][B->len_aln]='\0';
5280   return B;
5281 }
5282 Alignment * or_aln2pos_aln (Alignment *A, Alignment *S, char *mode)
5283 {
5284   Alignment *B;
5285
5286   int *pos;
5287   char outaln[100], outtree[100];
5288
5289
5290   strget_param (mode, "_OUTALN_", "NO", "%s", outaln);
5291   strget_param (mode, "_OUTTREE_", "NO", "%s", outtree);
5292
5293   set_sar (A, S, mode);
5294   pos=aln2predictive_positions (A, S,mode, PRINT);
5295
5296   B=display_pos (A, S, pos, mode);
5297
5298
5299   if (!strm(outaln, "NO")) vfclose (output_aln (B, vfopen (outaln, "w")));
5300   if (!strm(outtree, "NO"))vfclose (print_tree (aln2tree(B), "newick", vfopen (outtree, "w")));
5301
5302   return B;
5303 }
5304 Alignment * or_sim(Alignment *A, Alignment *S, char *mode)
5305 {
5306   //Predict all the sequences that are not both in inS and inA
5307   int *pos;
5308
5309   set_sar (A, S, mode);
5310   pos=aln2predictive_positions (A, S,mode, PRINT);
5311   fprintf ( stdout, "R: %.3f", pos2sim (A,S, pos));
5312
5313   exit (EXIT_SUCCESS);
5314   return A;
5315 }
5316 ORP* or_self_predict(Alignment *A, Alignment *S, char *mode,int *pos, int print)
5317 {
5318   //Predict all the sequences that are not both in inS and inA
5319   Alignment *P;
5320   char ***motifs;
5321
5322
5323   int a;
5324
5325   int pre_set_pos=0;
5326   ORP *PR;
5327
5328
5329   set_sar (A, S, mode);
5330   PR=declare_or_prediction (S->nseq, A->nseq, A->len_aln);
5331   sprintf (PR->mode, "self");
5332   PR->A=A;
5333   PR->S=S;
5334
5335   if (!pos)
5336     {
5337       pos=aln2predictive_positions (A, S,mode,print);
5338       pre_set_pos=0;
5339     }
5340   else
5341     pre_set_pos=1;
5342
5343   for (a=0; a< A->len_aln; a++)
5344     PR->pos[a]=pos[a];
5345
5346
5347   PR->motif=motifs=compounds2motifs (A, S, pos,0, mode, print);
5348   P=PR->P=aln2prediction (A, motifs, pos);
5349
5350   if (!pre_set_pos)vfree (pos);
5351
5352   PR=new_evaluate_prediction (PR, mode,print);
5353   return PR;
5354 }
5355
5356
5357 Alignment * or_predict(Alignment *inA, Alignment *inS, char *mode)
5358 {
5359   //Predict all the sequences that are not both in inS and inA
5360   Alignment *P, *A, *S, *T;
5361   char ***motifs;
5362   int *pos;
5363
5364   int a, b;
5365
5366
5367
5368
5369
5370
5371   A=copy_aln (inA, NULL);
5372   S=copy_aln (inS, NULL);
5373   set_sar (A, S, mode);
5374
5375   pos=aln2predictive_positions (A, S,mode,PRINT);
5376   motifs=compounds2motifs (A, S, pos,0, mode, PRINT);
5377   T=get_prediction_target (inA, inS, mode);
5378
5379
5380   P=aln2prediction (T, motifs, pos);
5381   //recall=evaluate_prediction (S, P, mode);
5382   for ( a=0; a<P->len_aln; a++)
5383     {
5384       for (b=0; b<P->nseq; b++)
5385         {
5386           if (tolower(P->seq_al[b][a])=='i')fprintf (stdout, "\n>%20s %20s %c", T->name [0],get_compound_name (a, mode), P->seq_al[b][a]);
5387         }
5388     }
5389   fprintf ( stdout, "\n");
5390   return P;
5391 }
5392
5393 Alignment *get_prediction_target (Alignment *A, Alignment *S, char *param)
5394 {
5395   char **name;
5396   int n, a;
5397   Alignment *T;
5398
5399   T=copy_aln (A, NULL);
5400   name=declare_char (A->nseq, 100);
5401   for (n=0,a=0; a< A->nseq; a++)
5402     {
5403       if ( name_is_in_list (A->name[a], S->name, S->nseq, 100)==-1)
5404         {
5405           sprintf (name[n++], "%s", A->name[a]);
5406         }
5407     }
5408   T=reorder_aln (T,name, n);
5409   return T;
5410 }
5411
5412 Alignment *set_sar (Alignment *A, Alignment *S, char *param)
5413 {
5414   char **name;
5415   int n, a;
5416
5417   name=declare_char (A->nseq, 100);
5418   for (n=0,a=0; a< A->nseq; a++)
5419     {
5420       if ( name_is_in_list (A->name[a], S->name, S->nseq, 100)!=-1)
5421         {
5422           sprintf (name[n++], "%s", A->name[a]);
5423         }
5424     }
5425   A=reorder_aln (A,name, n);
5426   S=reorder_aln (S,name, n);
5427   free_char (name, -1);
5428   return S;
5429 }
5430
5431 ORP* new_evaluate_prediction  (ORP *PR, char *mode, int print)
5432 {
5433   int a,b, i, r, p;
5434   int tp, tn, fp, fn;
5435   float sn, sp, sen2, best;
5436   float tot_best_seq=0;
5437   float tot_best_comp=0;
5438   Alignment *P, *R;
5439
5440   int ns=0;
5441   float *recall;
5442
5443
5444   P=PR->P;
5445   R=PR->S;
5446
5447   recall=vcalloc (P->len_aln, sizeof (float));
5448   if (P->len_aln!=R->len_aln)
5449     {
5450       HERE ("Mismatch between number of compounds in prediction and reference");
5451       exit (0);
5452     }
5453   if (print==PRINT)fprintf ( stdout, "\n");
5454
5455   for (a=0; a<P->nseq; a++)
5456     {
5457       tp=tn=fp=fn=0;
5458       if ((i=name_is_in_list (P->name[a], R->name, R->nseq, 100))!=-1)
5459         {
5460
5461           for (b=0;b<P->len_aln; b++)
5462             {
5463               r=R->seq_al[i][b];
5464               p=P->seq_al[a][b];
5465
5466               if ( p=='I' && r=='I')tp++;
5467               else if ( p=='I' && r=='O')fp++;
5468               else if ( p=='O' && r=='I')fn++;
5469               else if ( p=='O' && r=='O')tn++;
5470             }
5471           rates2sensitivity (tp, tn, fp, fn, &sp, &sn, &sen2, &best);
5472           if (print==PRINT)fprintf (stdout, ">%-s sp: %.2f sn: %.2f sn2: %.2f best: %.2f [SEQ]\n",P->name[a], sp, sn, sen2, best);
5473           if ( best>0)
5474             {
5475               ns++;
5476               tot_best_seq+=best;
5477             }
5478         }
5479     }
5480   if (ns)
5481     {
5482       tot_best_seq/=ns;
5483     }
5484   if (print==PRINT)fprintf ( stdout, ">TotSeq sp: %.2f N: %d[SEQ]\n",tot_best_seq, ns);
5485
5486   tot_best_comp=0;
5487   for (ns=0,b=0; b<P->len_aln; b++)
5488     {
5489       tp=tn=fp=fn=0;
5490       for (a=0; a<P->nseq;a++)
5491         {
5492           if ((i=name_is_in_list (P->name[a], R->name, R->nseq, 100))!=-1)
5493             {
5494               r=R->seq_al[i][b];
5495               p=P->seq_al[a][b];
5496
5497               if ( p=='I' && r=='I'){PR->tp++;tp++;}
5498               else if ( p=='I' && r=='O'){PR->fp++;fp++;}
5499               else if ( p=='O' && r=='I'){PR->fn++;fn++;}
5500               else if ( p=='O' && r=='O'){PR->tn++;tn++;}
5501             }
5502         }
5503       rates2sensitivity (tp, tn, fp, fn, &sp, &sn, &sen2, &best);
5504
5505       if (print==PRINT) fprintf (stdout, ">%-25s sp: %.2f sn: %.2f sen2: %.2f best: %.2f [COMP]\n",get_compound_name (b, mode), PR->sp, PR->sn, PR->sen2,PR->best);
5506       if ( best>0)
5507         {
5508           ns++;
5509           tot_best_comp+=best;
5510         }
5511     }
5512
5513   if (ns)
5514     {
5515       tot_best_comp/=ns;
5516     }
5517   rates2sensitivity (PR->tp, PR->tn, PR->fp,PR->fn,&(PR->sp), &(PR->sn), &(PR->sen2), &(PR->best));
5518   if (print==PRINT)fprintf ( stdout, ">FullTot sp: %.2f sn: %.2f sen2: %.2f best: %.2f  N: %d[COMP]\n", PR->sp, PR->sn, PR->sen2,PR->best, ns);
5519   P->score_aln=(int)((float)1000*(PR->best));
5520   return PR;
5521 }
5522 float  evaluate_prediction  (Alignment *R, Alignment *P, char *mode, int print)
5523 {
5524   int a,b, i, r, p;
5525   int tp, tn, fp, fn;
5526   int tot_tp, tot_tn, tot_fp, tot_fn;
5527   float sn, sp, sen2, best;
5528   float tot_sp=0;
5529   float tot_sn=0;
5530   float tot_sen2=0;
5531   float tot_best_seq=0;
5532   float tot_best_comp=0;
5533   float tot_best=0;
5534
5535   int ns=0;
5536   float *recall;
5537
5538
5539
5540
5541   recall=vcalloc (P->len_aln, sizeof (float));
5542   if (P->len_aln!=R->len_aln)
5543     {
5544       HERE ("Mismatch between number of compounds in prediction and reference");
5545       exit (0);
5546     }
5547   if (print==PRINT)fprintf ( stdout, "\n");
5548   for (a=0; a<P->nseq; a++)
5549     {
5550       tp=tn=fp=fn=0;
5551       if ((i=name_is_in_list (P->name[a], R->name, R->nseq, 100))!=-1)
5552         {
5553
5554           for (b=0;b<P->len_aln; b++)
5555             {
5556               r=R->seq_al[i][b];
5557               p=P->seq_al[a][b];
5558
5559               if ( p=='I' && r=='I')tp++;
5560               else if ( p=='I' && r=='O')fp++;
5561               else if ( p=='O' && r=='I')fn++;
5562               else if ( p=='O' && r=='O')tn++;
5563             }
5564           rates2sensitivity (tp, tn, fp, fn, &sp, &sn, &sen2, &best);
5565           if (print==PRINT)fprintf (stdout, ">%-s sp: %.2f sn: %.2f sn2: %.2f best: %.2f [SEQ]\n",P->name[a], sp, sn, sen2, best);
5566           if ( best>0)
5567             {
5568               ns++;
5569               tot_best_seq+=best;
5570               tot_sn+=sn;
5571               tot_sp+=sp;
5572               tot_sen2+=sen2;
5573             }
5574         }
5575     }
5576   if (ns)
5577     {
5578       tot_best_seq/=ns;
5579       tot_sn/=ns;
5580       tot_sp/=ns;
5581       tot_sen2/=ns;
5582     }
5583   if (print==PRINT)fprintf ( stdout, ">Tot sp: %.2f sn: %.2f sen2: %.2f best: %.2f  N: %d[SEQ]\n", tot_sp, tot_sn, tot_sen2,tot_best_seq, ns);
5584
5585   tot_fp=tot_fn=tot_tp=tot_tn=0;
5586   tot_sp=tot_sn=tot_sen2=tot_best_comp=0;
5587   for (ns=0,b=0; b<P->len_aln; b++)
5588     {
5589       tp=tn=fp=fn=0;
5590       for (a=0; a<P->nseq;a++)
5591         {
5592           if ((i=name_is_in_list (P->name[a], R->name, R->nseq, 100))!=-1)
5593             {
5594               r=R->seq_al[i][b];
5595               p=P->seq_al[a][b];
5596
5597               if ( p=='I' && r=='I'){tot_tp++;tp++;}
5598               else if ( p=='I' && r=='O'){tot_fp++;fp++;}
5599               else if ( p=='O' && r=='I'){tot_fn++;fn++;}
5600               else if ( p=='O' && r=='O'){tot_tn++;tn++;}
5601             }
5602         }
5603       rates2sensitivity (tp, tn, fp, fn, &sp, &sn, &sen2, &best);
5604
5605       if (print==PRINT) fprintf (stdout, ">%-25s sp: %.2f sn: %.2f sen2: %.2f best: %.2f [COMP]\n",get_compound_name (b, mode), sp, sn, sen2,best);
5606       recall[b]=sen2;
5607       if ( best>0)
5608         {
5609           ns++;
5610           tot_best_comp+=best;
5611           tot_sn+=sn;
5612           tot_sp+=sp;
5613           tot_sen2+=sen2;
5614         }
5615     }
5616
5617   if (ns)
5618     {
5619       tot_best_comp/=ns;
5620       tot_sn/=ns;
5621       tot_sp/=ns;
5622       tot_sen2/=ns;
5623     }
5624   rates2sensitivity (tot_tp, tot_tn, tot_fp,tot_fn,&tot_sp, &tot_sn, &tot_sen2, &tot_best);
5625   if (print==PRINT)fprintf ( stdout, ">FullTot sp: %.2f sn: %.2f sen2: %.2f best: %.2f  N: %d[COMP]\n", tot_sp, tot_sn, tot_sen2,tot_best, ns);
5626   return tot_best;
5627 }
5628
5629
5630
5631 Alignment * aln2prediction (Alignment *A,char ***motif, int *pos)
5632 {
5633   int a, b,nc, nl;
5634   int *list;
5635   char **array, **sar;
5636   Alignment *R;
5637   Sequence *S;
5638   nc=read_array_size ((void *)motif, sizeof (char***));
5639
5640
5641   list=pos2list (pos, A->len_aln, &nl);
5642
5643
5644   array=declare_char (A->nseq, nl+1);
5645   sar=declare_char(A->nseq, nc+1);
5646   for (a=0; a<A->nseq; a++)
5647     {
5648       for (b=0; b<nl; b++)
5649         array[a][b]=A->seq_al[a][list[b]];
5650     }
5651
5652   for (a=0; a<nc; a++)
5653     {
5654       for (b=0; b<A->nseq; b++)
5655         {
5656
5657           sar[b][a]=(match_motif (array[b], motif[a]))?'I':'O';
5658         }
5659     }
5660
5661
5662   S=fill_sequence_struc (A->nseq,sar,A->name);
5663   R=seq2aln (S, NULL, KEEP_GAP);
5664   free_sequence (S, S->nseq);
5665   free_char (sar, -1);
5666   vfree (list);
5667   free_char (array, -1);
5668   return R;
5669 }
5670
5671 int *   file2pos_list (Alignment *A, char *posfile)
5672 {
5673   char ***file;
5674   int **index;
5675   int *pos;
5676   int i, n, p;
5677
5678   //pos_file:<P: seqname pos score\n>
5679   //          1  2       3   4
5680
5681
5682   if ( !check_file_exists (posfile))
5683     {
5684       printf_exit ( EXIT_FAILURE, stderr, "ERROR: Could not read posfile %s\n", posfile);
5685     }
5686
5687   file=file2list (posfile, " ");
5688
5689   index=aln2inv_pos (A);
5690   pos=vcalloc ( A->len_aln, sizeof (int));
5691
5692   n=0;
5693   while (file[n])
5694     {
5695
5696       if ( !strm (file[n][1], "P:"));
5697       else
5698         {
5699           if ( (strm (file[n][2], "cons")))
5700             p=atoi(file[n][3])-1;
5701           else
5702             {
5703               i=name_is_in_list ( file[n][2], A->name, A->nseq, MAXNAMES+1);
5704               if (i!=-1)
5705                 p=index[i][atoi(file[n][3])]-1;
5706               else p=-1;
5707             }
5708           if (p!=-1)pos[p]+=atoi(file[n][4]);
5709         }
5710       n++;
5711     }
5712
5713
5714   free_int (index, -1);
5715   free_arrayN ( (char **)file, 3);
5716   return pos;
5717 }
5718 int *   aln2predictive_positions (Alignment *A, Alignment *B, char *mode, int print)
5719 {
5720   char posmode[100];
5721
5722   if (!mode) return NULL;
5723   HERE ("%s", mode);
5724   strget_param (mode, "_POSMODE_", "scan", "%s", posmode);
5725   if ( strm (posmode, "mat"))return aln2predictive_positions_mat (A, B, mode, print);
5726   else if ( strm (posmode, "scan")) return aln2predictive_positions_scan (A, B, mode, print);
5727   else
5728     {
5729       printf_exit (EXIT_FAILURE,stderr, "ERROR: %s is an unknown _POSMODE_ mode",posmode);
5730       return NULL;
5731     }
5732 }
5733
5734 int *   aln2predictive_positions_mat (Alignment *A, Alignment *B, char *mode, int print)
5735   {
5736     int a, b, c,gap,  res1, res2, sar1, sar2, npos, s, idscore;
5737     float id1,id2,id3,nid1,nid2,nid3;
5738     int **pos, *fpos;
5739     pos=declare_int (A->len_aln,2);
5740     fpos=vcalloc ( A->len_aln, sizeof (int));
5741
5742     strget_param (mode, "_NPOS_", "2", "%d", &npos);
5743     for ( a=0; a< A->len_aln; a++)
5744       {
5745         pos[a][0]=a;
5746         id1=id2=id3=nid1=nid2=nid3=0;
5747         for ( gap=0,b=0; b<A->nseq; b++)gap+=(A->seq_al[b][a]=='-');
5748         if ( gap>0){pos[a][1]=0;continue;}
5749
5750         for (s=0; s<B->len_aln; s++)
5751           {
5752             for ( gap=0,b=0; b<A->nseq-1; b++)
5753               {
5754                 sar1=B->seq_al[b][s];
5755                 res1=A->seq_al[b][a];
5756
5757                 for ( c=b+1; c<A->nseq; c++)
5758                   {
5759                     sar2=B->seq_al[c][s];
5760                     res2=A->seq_al[c][a];
5761
5762                     idscore=(res1==res2)?1:0;
5763                     if ( sar1 == 'I' && sar2=='I'){id1+=idscore;nid1++;}
5764                     else if ( sar1 =='0' && sar2=='0'){id2+=idscore;nid2++;}
5765                     else {id3+=idscore; nid3++;}
5766
5767                   }
5768               }
5769             id1=(nid1==0)?1:id1/nid1;
5770             id2=(nid1==0)?1:id2/nid2;
5771             id3=(nid3==0)?1:id3/nid3;
5772             pos[a][1]=(int)((float)1000*id1*(1-id3));
5773
5774           }
5775       }
5776
5777     sort_int (pos, 2,1, 0, A->len_aln-1);
5778     for ( a=MAX(0,(A->len_aln-npos));a<A->len_aln; a++)
5779       {
5780         fpos[pos[a][0]]=1;
5781       }
5782
5783     free_int (pos, -1);
5784     return fpos;
5785   }
5786 int *   aln2predictive_positions_scan (Alignment *A, Alignment *B, char *mode, int print)
5787 {
5788   int a, b, c, best_pos,nl, nplist=0, max, posw;
5789   float best_score, score;
5790   static int *list, *tpos,**plist,*array;
5791   int *pos;
5792
5793
5794   char posfile[100];
5795   char predmode[100];
5796   char target_posfile[100];
5797
5798
5799
5800   if (!A)
5801     {
5802       vfree (list);
5803       vfree (tpos);
5804
5805       free_int (plist, -1);
5806       vfree (array);
5807       return NULL;
5808     }
5809
5810   strget_param (mode, "_PREDMODE_", "ID", "%s", predmode);
5811   strget_param (mode, "_POSW_", "1", "%d", &posw);
5812   strget_param (mode, "_NPOS_", "2", "%d", &max);
5813   strget_param (mode, "_POSFILE_", "NO", "%s", posfile);
5814   strget_param (mode, "_TPOSFILE_", "NO", "%s", target_posfile);
5815
5816   if ( !strm(posfile, "NO"))return file2pos_list (A,posfile);
5817   if ( !strm(target_posfile, "NO"))tpos=file2pos_list (A,target_posfile);
5818   else
5819     {
5820       tpos=vcalloc (A->len_aln, sizeof (int));
5821       for (a=0; a<A->len_aln; a++)tpos[a]=1;
5822     }
5823
5824   //Declare the positions that are going to be scanned
5825
5826
5827   if (posw==1)
5828     {
5829       plist=declare_int (A->len_aln, 2);
5830       nplist=0;
5831       for (a=0; a<A->len_aln; a++)
5832         {
5833           if(tpos[a])
5834             {
5835               plist[nplist][0]=1;
5836               plist[nplist][1]=a;
5837               nplist++;
5838             }
5839         }
5840     }
5841   else if ( posw==2)
5842     {
5843       nplist=0;
5844       plist=declare_int (A->len_aln*A->len_aln, 3);
5845       for (a=0; a<A->len_aln; a++)
5846         for (b=0; b<A->len_aln; b++)
5847           {
5848             plist[nplist][1]=a;
5849             plist[nplist][2]=b;
5850             plist[nplist][0]=2;
5851             nplist++;
5852           }
5853     }
5854   else if ( posw==3)
5855     {
5856       nplist=0;
5857       plist=declare_int (A->len_aln*A->len_aln*A->len_aln, 3);
5858       for (a=0; a<A->len_aln; a++)
5859         for (b=0; b<A->len_aln; b++)
5860           {
5861             plist[nplist][1]=a;
5862             plist[nplist][2]=b;
5863             plist[nplist][3]=c;
5864
5865
5866             plist[nplist][0]=3;
5867             nplist++;
5868           }
5869     }
5870
5871
5872   pos=vcalloc ( A->len_aln, sizeof (int));
5873   if (max==0)max=A->len_aln;
5874   else if ( max==-1)
5875     {
5876       for (a=0; a<A->len_aln; a++)if (tpos[a]){pos[a]=1;}
5877       aln2predictive_positions_scan (NULL, NULL, NULL, 0);
5878       return pos;
5879     }
5880
5881
5882
5883   pos=vcalloc ( A->len_aln, sizeof (int));
5884   list=vcalloc (A->len_aln, sizeof (int));
5885   nl=0;
5886
5887
5888
5889   for (a=0; a< max; a++)
5890     {
5891       int previous_best_pos=-1;
5892       for (best_score=-9999,best_pos=0,b=0; b<nplist; b++)
5893         {
5894           for (c=0; c<nl; c++)pos[list[c]]=1;
5895           for (c=1; c<=plist[b][0]; c++)pos[plist[b][c]]=1;
5896
5897           if (strm(predmode, "R"))score=sar_aln2r(A,B,pos,0);
5898           else if ( strm (predmode, "ID"))
5899             {
5900               score  =or_id_evaluate (A, B, mode, pos,NO_PRINT);
5901             }
5902           else if ( strm (predmode, "BP2"))score =or_loo_evaluate2 (A, B, mode, pos,NO_PRINT);
5903           else
5904             {
5905               HERE ("Unknown mode: %s", predmode);
5906               exit (0);
5907             }
5908           if ( score>best_score)
5909             {
5910               best_score=score;
5911               best_pos=b;
5912             }
5913           for (c=1; c<=plist[b][0]; c++)pos[plist[b][c]]=0;
5914
5915         }
5916       if (best_pos==previous_best_pos)break;
5917       else previous_best_pos=best_pos;
5918
5919       //update the best_pos_list
5920       for (b=1; b<=plist[best_pos][0]; b++)
5921         list[nl++]=plist[best_pos][b];
5922
5923
5924       if ( print==PRINT)
5925         {
5926           for (b=0; b<nl; b++) pos[list[b]]=1;
5927           fprintf ( stdout, "\nP_current: ");
5928           for ( c=1; c<=plist[best_pos][0]; c++)fprintf ( stdout, "%d ",plist[best_pos][c]+1);
5929           fprintf ( stdout, " S: %.3f D: %d R:%d", best_score, (int)sar_aln2delta(A,B, pos,0), nl);
5930         }
5931       for (b=0; b<nl; b++) pos[list[b]]=0;
5932     }
5933
5934   for (a=0; a<nl; a++)
5935     pos[list[a]]=1;
5936   if (print==PRINT)fprintf ( stdout, "\nR_best: %.3f with %d pos" ,best_score, nl);
5937
5938   aln2predictive_positions_scan (NULL, NULL, NULL, 0);
5939
5940   return pos;
5941 }
5942
5943 char *** compounds2motifs (Alignment *A, Alignment *B, int *pos, int depth, char *mode, int print)
5944 {
5945   char ***motifs;
5946   int a;
5947
5948   motifs=vcalloc (B->len_aln, sizeof (char**));
5949   for (a=0; a<B->len_aln; a++)
5950     {
5951
5952       motifs[a]=compound2motif (A, B, pos, depth, a, mode, print);
5953     }
5954
5955   return motifs;
5956 }
5957 char ** compound2regexp_motif (Alignment *A, Alignment *B, int *pos, int depth, int c, char *mode, int print);
5958 char ** compound2word_motif (Alignment *A, Alignment *B, int *pos, int depth, int c, char *mode, int print);
5959
5960 char ** compound2motif (Alignment *A, Alignment *B, int *pos, int depth, int c, char *mode, int print)
5961 {
5962   char mmode[100];
5963
5964   strget_param (mode, "_MOTIFMODE_", "word", "%s", mmode); //words, regexp
5965   if ( strm (mmode, "regexp"))return compound2regexp_motif (A,B,pos, depth, c, mode, print);
5966   else if ( strm (mmode, "word"))return compound2word_motif (A,B,pos, depth, c, mode, print);
5967  else return NULL;}
5968 char ** compound2word_motif (Alignment *A, Alignment *B, int *pos, int depth, int c, char *mode, int print)
5969 {
5970   int a,l;
5971   char *word, **motif;
5972   float score;
5973
5974
5975   word=or_id_evaluate2 (A, B, mode, pos,print, &score);
5976   if ( !word) return NULL;
5977   l=strlen (word);
5978
5979   motif=declare_char (l+1, 2);
5980   for (a=0; a<l; a++)motif[a][0]=word[a];
5981
5982   if (print==PRINT)
5983     {
5984       fprintf ( stdout, "\nMotifCompound %25s best: %.2f  motif: ", get_compound_name(c, mode),score);
5985       for (a=0; a<l; a++)
5986         {
5987           fprintf ( stdout, "[%2s]",motif[a]);
5988         }
5989     }
5990   vfree (word);
5991   return motif;
5992 }
5993
5994
5995
5996 char ** compound2regexp_motif (Alignment *A, Alignment *B, int *pos, int depth, int c, char *mode, int print)
5997 {
5998   static Alignment *I;
5999   static Alignment *O;
6000   int a, b, o, i;
6001
6002   float tp,tn,fp,fn,best, sp, sn, sen2, best_sn, best_sp, best_sen2;
6003   float best_pred=-1;
6004   int best_motif=0;
6005
6006
6007
6008   char ***alp=NULL;
6009   int *alp_size=NULL;
6010
6011   char *motif_file;
6012   int n;
6013   char **m, **m2;
6014
6015   char *buf=NULL;
6016   char *best_buf;
6017   FILE *fpp;
6018
6019   free_aln (I);
6020   free_aln (O);
6021
6022   I=copy_aln(A, NULL);
6023   O=copy_aln(A, NULL);
6024
6025   if (depth==0)
6026     strget_param (mode, "_DEPTH_", "2", "%d", &depth);
6027
6028   I->nseq=O->nseq=I->len_aln=O->len_aln=0;
6029   for (a=0; a<A->len_aln; a++)
6030     {
6031       if (pos[a])
6032         {
6033           for (i=o=0,b=0; b<A->nseq; b++)
6034             {
6035               if ( is_gap(A->seq_al[b][a]))return 0;
6036               if (B->seq_al[b][c]=='I')I->seq_al[i++][I->len_aln]=A->seq_al[b][a];
6037               else O->seq_al[o++][O->len_aln]=A->seq_al[b][a];
6038             }
6039           I->len_aln++;
6040           O->len_aln++;
6041         }
6042     }
6043
6044   if (O->len_aln==0 || I->len_aln==0) return 0;
6045   O->nseq=o;
6046   I->nseq=i;
6047   for (a=0; a<o; a++)O->seq_al[a][O->len_aln]='\0';
6048   for (a=0; a<i; a++)I->seq_al[a][I->len_aln]='\0';
6049
6050   if (!I->nseq) return NULL;
6051
6052
6053
6054   best_pred=best_motif=best_sn=best_sp=best_sen2=0;
6055
6056   motif_file=vtmpnam (NULL);
6057
6058   n=0;
6059   if (depth>0)
6060     {
6061        alp=vcalloc ( sizeof (char**), I->len_aln);
6062        alp_size= vcalloc ( I->len_aln, sizeof (int));
6063        for (a=0; a<I->len_aln; a++)
6064          {
6065            char *col;
6066            alp[a]=string2alphabet ( (col=aln_column2string (I,a)),depth, &alp_size[a]);
6067            vfree (col);
6068          }
6069        generate_array_string_list (I->len_aln, alp, alp_size, &n, motif_file, OVERLAP);
6070     }
6071   else
6072     {
6073       int *used;
6074       char r;
6075
6076       used=vcalloc (256, sizeof (int));
6077       fpp=vfopen (motif_file,"w");
6078       for (a=0;a<I->len_aln; a++)
6079         {
6080           for (b=0; b<I->nseq; b++)
6081             {
6082               r=I->seq_al[b][a];
6083               if (!used[(int)r]){fprintf (fpp, "%c", r);used[(int)r]=1;}
6084             }
6085           for (b=0; b<I->nseq; b++)
6086             {
6087               r=I->seq_al[b][a];
6088               used[(int)r]=0;
6089             }
6090           fprintf (fpp, " ");
6091         }
6092       fprintf (fpp, "\n");
6093       vfree (used);
6094       vfclose (fpp);
6095
6096       n=1;
6097       depth=I->nseq;
6098     }
6099
6100   buf=vcalloc (2*(I->len_aln*depth)+1, sizeof (char));
6101   best_buf=vcalloc (2*(I->len_aln*depth)+1, sizeof (char));
6102   fpp=vfopen (motif_file, "r");
6103
6104   for (a=0; a<n; a++)
6105     {
6106
6107       buf=vfgets (buf, fpp);
6108       m2=string2list (buf);
6109       m2++;
6110
6111
6112       tp=tn=fp=fn=0;
6113
6114       for (b=0; b<I->nseq; b++)
6115         {
6116           if (match_motif (I->seq_al[b], m2))tp++;
6117           else fn++;
6118         }
6119       for (b=0; b<O->nseq; b++)
6120         {
6121           if (match_motif (O->seq_al[b], m2))fp++;
6122           else tn++;
6123         }
6124       rates2sensitivity (tp, tn, fp, fn, &sp, &sn, &sen2, &best);
6125
6126       if (best>= best_pred)
6127         {
6128           best_pred=best;
6129           best_sp=sp;
6130           best_sen2=sen2;
6131           best_sn=sn;
6132           sprintf (best_buf, "%s", buf);
6133         }
6134       m2--;
6135       free_char (m2, -1);
6136     }
6137   vfclose (fpp);
6138   if (print==PRINT)fprintf ( stdout, "\nMotifCompound %25s sp: %.2f sn: %.2f sen2: %.2f best: %.2f  motif: ", get_compound_name(c, mode), best_sp, best_sn, best_sen2, best_pred);
6139   m2=string2list (best_buf);
6140   m=declare_char (I->len_aln+1, depth+1);
6141
6142   for (a=0; a<I->len_aln; a++)
6143         {
6144           sprintf (m[a], "%s", m2[a+1]);
6145           if (print==PRINT) fprintf ( stdout, "[%2s]",m[a]);
6146         }
6147   if (print==PRINT)fprintf ( stdout, " N-motifs %d", n);
6148   free_char (m2, -1);
6149
6150   if (alp)free_arrayN((void ***) alp, 3);
6151   if (alp_size)vfree (alp_size);
6152   vfree (buf); vfree(best_buf);
6153
6154   return m;
6155 }
6156
6157 double pos2sim (Alignment *A, Alignment *B, int *pos)
6158 {
6159   return sar_aln2r (A, B,pos, PRINT);
6160 }
6161 double  sar_aln2r (Alignment *A, Alignment *B, int *pos, int print)
6162 {
6163   int a, b, c, d,r1, r2, n, score, sim;
6164   double *r, result;
6165   static double **slist;
6166   int declare=0;
6167   static int **M;
6168
6169
6170
6171   if (!M)M=read_matrice ("blosum62mt");
6172   if (!slist)
6173     {
6174       int maxslist;
6175       maxslist=A->nseq*A->nseq*10;
6176       slist=declare_double (maxslist, 2);
6177     }
6178
6179   if (pos==NULL)
6180     {
6181
6182       declare=1;
6183       pos=vcalloc ( A->len_aln+1, sizeof (int));
6184       for (a=0; a<A->len_aln; a++)pos[a]=1;
6185       pos[a]=-1;
6186
6187     }
6188
6189   for (n=0,a=0; a< A->nseq-1; a++)
6190     {
6191
6192       for (b=a+1; b<A->nseq; b++)
6193         {
6194
6195
6196           for (sim=d=0,c=0; c<A->len_aln; c++)
6197             {
6198
6199               if (pos[c]==0)continue;
6200
6201               r1=A->seq_al[a][c];
6202               r2=A->seq_al[b][c];
6203               if (is_gap(r1) || is_gap(r2))return 0;
6204
6205               sim+=M[r1-'A'][r2-'A']*pos[c];
6206               d+=MAX((M[r1-'A'][r1-'A']),(M[r2-'A'][r2-'A']));
6207             }
6208           sim=(d==0)?0:(100*sim)/d;
6209           score=(int)get_sar_sim(B->seq_al[a], B->seq_al[b]);
6210           slist[n][0]=(double)sim;
6211           slist[n][1]=(double)score;
6212           if (print==PRINT)fprintf ( stdout, "SIM: %d %d [%s %s]\n", sim, score, A->name[a], A->name[b]);
6213           n++;
6214         }
6215     }
6216
6217   r=return_r(slist, n);
6218   for (a=0; a<n; a++)slist[a][0]=slist[a][1]=0;
6219   result=r[0];
6220   vfree (r);
6221   if (declare) vfree (pos);
6222
6223   return result;
6224 }
6225
6226
6227 double sar_aln2delta (Alignment *A, Alignment *B, int *pos, int print)
6228 {
6229   static Alignment *I;
6230   static Alignment *O;
6231   int a, b, c, o, i;
6232   double delta=0;
6233   if (!I)
6234     {
6235       I=copy_aln(A, NULL);
6236       O=copy_aln(A, NULL);
6237     }
6238
6239
6240
6241   for (c=0; c<B->len_aln; c++)
6242     {
6243
6244       I->nseq=O->nseq=I->len_aln=O->len_aln=0;
6245       for (a=0; a<A->len_aln; a++)
6246         {
6247           if (pos[a])
6248             {
6249               for (i=o=0,b=0; b<B->nseq; b++)
6250                 {
6251                   if ( is_gap(A->seq_al[b][a]))return 0;
6252                   if (B->seq_al[b][c]=='I')I->seq_al[i++][I->len_aln]=A->seq_al[b][a];
6253                   else O->seq_al[o++][O->len_aln]=A->seq_al[b][a];
6254                 }
6255               I->len_aln++;
6256               O->len_aln++;
6257             }
6258         }
6259       if (O->len_aln==0 || I->len_aln==0) return 0;
6260       O->nseq=o;
6261       I->nseq=i;
6262       for (a=0; a<o; a++)O->seq_al[a][O->len_aln]='\0';
6263       for (a=0; a<i; a++)I->seq_al[a][I->len_aln]='\0';
6264
6265       delta+=aln2sim(I,"blosum62mt")-aln2sim(O, "blosum62mt");
6266
6267     }
6268
6269   return delta;
6270 }
6271
6272 char * get_compound_name (int c, char *mode)
6273 {
6274   static int isset;
6275   static Alignment *S;
6276   static char *lname;
6277
6278   if (!isset)
6279     {
6280       char *comp_list;
6281       isset=1;
6282       lname=vcalloc (100, sizeof (char));
6283
6284       if (!mode);
6285       else
6286         {
6287           strget_param (mode, "_COMPLIST_", "NO", "%s", comp_list=vcalloc (100, sizeof (char)));
6288           if (strm(comp_list, "NO"));
6289           else
6290             {
6291               S=main_read_aln (comp_list, NULL);
6292               vfree (comp_list);
6293             }
6294         }
6295     }
6296   if (!S || c>=S->nseq)sprintf (lname, "%d", c);
6297   else
6298     {
6299       sprintf (lname, "%s", S->name [c]);
6300     }
6301   return lname;
6302 }
6303 ORP * declare_or_prediction ( int ncomp, int nseq, int len)
6304 {
6305   ORP *P;
6306   P=vcalloc ( 1, sizeof (ORP));
6307   P->ncomp=ncomp;
6308   P->nseq=nseq;
6309   P->len=len;
6310   P->PR=NULL;
6311
6312   P->pos=vcalloc (len+1, sizeof (int));
6313
6314   return P;
6315 }
6316
6317 void free_orp_list ( ORP**P)
6318 {
6319   int a=0;
6320   while (P[a])
6321     {
6322       free_orp(P[a++]);
6323     }
6324 }
6325 void free_orp ( ORP*P)
6326 {
6327   if (!P) return;
6328   free_aln (P->A);
6329   free_aln (P->S);
6330   free_aln (P->P);
6331   vfree (P->pos);
6332   free_arrayN((void **)P->motif, 3);
6333   if (P->PR)free_orp(P->PR);
6334   vfree (P);
6335 }
6336
6337
6338
6339
6340
6341
6342
6343
6344
6345
6346
6347
6348
6349
6350
6351 /*********************************COPYRIGHT NOTICE**********************************/
6352 /*© Centro de Regulacio Genomica */
6353 /*and */
6354 /*Cedric Notredame */
6355 /*Tue Oct 27 10:12:26 WEST 2009. */
6356 /*All rights reserved.*/
6357 /*This file is part of T-COFFEE.*/
6358 /**/
6359 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
6360 /*    it under the terms of the GNU General Public License as published by*/
6361 /*    the Free Software Foundation; either version 2 of the License, or*/
6362 /*    (at your option) any later version.*/
6363 /**/
6364 /*    T-COFFEE is distributed in the hope that it will be useful,*/
6365 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
6366 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
6367 /*    GNU General Public License for more details.*/
6368 /**/
6369 /*    You should have received a copy of the GNU General Public License*/
6370 /*    along with Foobar; if not, write to the Free Software*/
6371 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
6372 /*...............................................                                                                                      |*/
6373 /*  If you need some more information*/
6374 /*  cedric.notredame@europe.com*/
6375 /*...............................................                                                                                                                                     |*/
6376 /**/
6377 /**/
6378 /*      */
6379 /*********************************COPYRIGHT NOTICE**********************************/