binaries/src/tcoffee/t_coffee_source/evaluate_for_struc.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <stdarg.h>
   5 #include <ctype.h>
   6 #include <string.h>
   7
   8 #include "io_lib_header.h"
   9 #include "util_lib_header.h"
  10 #include "define_header.h"
  11 #include "dp_lib_header.h"
  12 void print_atom ( Atom*A);
  13
  14 float **** quantile_apdb_filtration ( Alignment *A, float ****residues, Constraint_list *CL,Pdb_param *PP, FILE *fp);
  15 float **** irmsdmin_apdb_filtration ( Alignment *A, float ****residues, Constraint_list *CL,Pdb_param *PP, FILE *fp);
  16 int apdb ( int argc, char *argv[])
  17     {
  18
  19         Constraint_list *CL=NULL;
  20         Sequence  *S=NULL;
  21         Alignment *A=NULL;
  22         Alignment *EA=NULL;
  23         Pdb_param *pdb_param;
  24
  25         Fname *F=NULL;
  26         char *file_name;
  27         int a,c;
  28
  29         int n_pdb;
  30
  31 /*PARAMETERS VARIABLES*/
  32         int garbage;
  33         char *parameters;
  34         FILE *fp_parameters;
  35
  36         int quiet;
  37         char *se_name;
  38         FILE *le=NULL;
  39
  40         char **list_file;
  41         int n_list;
  42         char **struc_to_use;
  43         int n_struc_to_use;
  44
  45         char *aln;
  46         char *repeat_seq;
  47         char *repeat_pdb;
  48
  49         char *color_mode;
  50         char *comparison_io;
  51
  52         int  n_excluded_nb;
  53
  54         float maximum_distance;
  55         float similarity_threshold;
  56         float   md_threshold;
  57
  58
  59         int print_rapdb;
  60
  61         char  *outfile;
  62         char  *run_name;
  63
  64         char  *apdb_outfile;
  65         char *cache;
  66
  67         char **out_aln_format;
  68         int    n_out_aln_format;
  69
  70         char *output_res_num;
  71         char *local_mode;
  72         float filter;
  73         int filter_aln;
  74         int irmsd_graph;
  75         int nirmsd_graph;
  76         int n_template_file;
  77         char **template_file_list;
  78         char *mode;
  79                 int prot_min_sim;
  80         int prot_max_sim;
  81         int prot_min_cov;
  82         int pdb_min_sim;
  83         int pdb_max_sim;
  84         int pdb_min_cov;
  85         int gapped;
  86
  87
  88         char *prot_blast_server;
  89         char *pdb_blast_server;
  90
  91
  92         char *pdb_db;
  93         char *prot_db;
  94         int min_ncol;
  95
  96         argv=standard_initialisation (argv, &argc);
  97
  98 /*PARAMETER PROTOTYPE:    READ PARAMETER FILE     */
  99                declare_name (parameters);
 100                get_cl_param(\
 101                             /*argc*/      argc             ,\
 102                             /*argv*/      argv             ,\
 103                             /*output*/    &le              ,\
 104                             /*Name*/      "-parameters"        ,\
 105                             /*Flag*/      &garbage        ,\
 106                             /*TYPE*/      "R_F"              ,\
 107                             /*OPTIONAL?*/ OPTIONAL         ,\
 108                             /*MAX Nval*/  1                ,\
 109                             /*DOC*/       "Read the files in the parameter file" ,\
 110                             /*Parameter*/ &parameters          ,\
 111                             /*Def 1*/     "NULL"             ,\
 112                             /*Def 2*/     "stdin"             ,\
 113                             /*Min_value*/ "any"            ,\
 114                             /*Max Value*/ "any"             \
 115                    );
 116                if ( parameters && parameters[0])
 117                   {
 118                   argv[argc]=vcalloc ( VERY_LONG_STRING, sizeof(char));
 119                   a=0;
 120                   fp_parameters=vfopen (parameters, "r");
 121                   while ((c=fgetc (fp_parameters))!=EOF)argv[1][a++]=c;
 122                   vfclose (fp_parameters);
 123                   argv[argc][a]='\0';
 124                   argc++;
 125                   argv=break_list ( argv, &argc, "=:;, \n");
 126                   }
 127 /*PARAMETER PROTOTYPE*/
 128      declare_name (se_name);
 129      get_cl_param(\
 130                             /*argc*/      argc          ,\
 131                             /*argv*/      argv          ,\
 132                             /*output*/    &le           ,\
 133                             /*Name*/      "-quiet"      ,\
 134                             /*Flag*/      &quiet        ,\
 135                             /*TYPE*/      "W_F"         ,\
 136                             /*OPTIONAL?*/ OPTIONAL      ,\
 137                             /*MAX Nval*/  1             ,\
 138                             /*DOC*/       "ND"          ,\
 139                             /*Parameter*/ &se_name      ,\
 140                             /*Def 1*/     "stderr"   ,\
 141                             /*Def 2*/     "/dev/null"   ,\
 142                             /*Min_value*/ "any"         ,\
 143                             /*Max Value*/ "any"          \
 144                    );
 145
 146      le=vfopen ( se_name, "w");
 147      fprintf ( le, "\nPROGRAM: %s\n",argv[0]);
 148
 149 /*PARAMETER PROTOTYPE:        IN */
 150         list_file=declare_char ( 200, STRING);
 151         n_list=get_cl_param(\
 152                             /*argc*/      argc          ,\
 153                             /*argv*/      argv          ,\
 154                             /*output*/    &le           ,\
 155                             /*Name*/      "-in"         ,\
 156                             /*Flag*/      &garbage      ,\
 157                             /*TYPE*/      "S"           ,\
 158                             /*OPTIONAL?*/ OPTIONAL      ,\
 159                             /*MAX Nval*/  200           ,\
 160                             /*DOC*/       "ND"          ,\
 161                             /*Parameter*/ list_file     ,\
 162                             /*Def 1*/    "",\
 163                             /*Def 2*/     "stdin"       ,\
 164                             /*Min_value*/ "any"         ,\
 165                             /*Max Value*/ "any"          \
 166                    );
 167 /*PARAMETER PROTOTYPE:        IN */
 168         struc_to_use=declare_char ( 200, STRING);
 169         n_struc_to_use=get_cl_param(\
 170                             /*argc*/      argc          ,\
 171                             /*argv*/      argv          ,\
 172                             /*output*/    &le           ,\
 173                             /*Name*/      "-struc_to_use"         ,\
 174                             /*Flag*/      &garbage      ,\
 175                             /*TYPE*/      "S"           ,\
 176                             /*OPTIONAL?*/ OPTIONAL      ,\
 177                             /*MAX Nval*/  200           ,\
 178                             /*DOC*/       "ND"          ,\
 179                             /*Parameter*/ struc_to_use     ,\
 180                             /*Def 1*/    "",\
 181                             /*Def 2*/     "stdin"       ,\
 182                             /*Min_value*/ "any"         ,\
 183                             /*Max Value*/ "any"          \
 184                    );
 185
 186 /*PARAMETER PROTOTYPE:        COMPARISON IO */
 187         declare_name (comparison_io);
 188         get_cl_param(\
 189                             /*argc*/      argc          ,\
 190                             /*argv*/      argv          ,\
 191                             /*output*/    &le           ,\
 192                             /*Name*/      "-io_format"  ,\
 193                             /*Flag*/      &garbage      ,\
 194                             /*TYPE*/      "S"           ,\
 195                             /*OPTIONAL?*/ OPTIONAL      ,\
 196                             /*MAX Nval*/  200           ,\
 197                             /*DOC*/       "ND"          ,\
 198                             /*Parameter*/  &comparison_io,\
 199                             /*Def 1*/    "hsgd0123456",\
 200                             /*Def 2*/     ""       ,\
 201                             /*Min_value*/ "any"         ,\
 202                             /*Max Value*/ "any"          \
 203                    );
 204 /*PARAMETER PROTOTYPE:        ALN */
 205         declare_name (aln);
 206         get_cl_param(\
 207                             /*argc*/      argc          ,\
 208                             /*argv*/      argv          ,\
 209                             /*output*/    &le           ,\
 210                             /*Name*/      "-aln"         ,\
 211                             /*Flag*/      &garbage      ,\
 212                             /*TYPE*/      "S"           ,\
 213                             /*OPTIONAL?*/ OPTIONAL      ,\
 214                             /*MAX Nval*/  1           ,\
 215                             /*DOC*/       "ND"          ,\
 216                             /*Parameter*/  &aln,\
 217                             /*Def 1*/    "",\
 218                             /*Def 2*/     "stdin"       ,\
 219                             /*Min_value*/ "any"         ,\
 220                             /*Max Value*/ "any"          \
 221                    );
 222 /*PARAMETER PROTOTYPE:        ALN */
 223
 224         declare_name (repeat_seq);
 225         get_cl_param(\
 226                             /*argc*/      argc          ,\
 227                             /*argv*/      argv          ,\
 228                             /*output*/    &le           ,\
 229                             /*Name*/      "-repeat_seq"         ,\
 230                             /*Flag*/      &garbage      ,\
 231                             /*TYPE*/      "S"           ,\
 232                             /*OPTIONAL?*/ OPTIONAL      ,\
 233                             /*MAX Nval*/  1           ,\
 234                             /*DOC*/       "ND"          ,\
 235                             /*Parameter*/  &repeat_seq,\
 236                             /*Def 1*/    "",\
 237                             /*Def 2*/     "stdin"       ,\
 238                             /*Min_value*/ "any"         ,\
 239                             /*Max Value*/ "any"          \
 240                    );
 241
 242 /*PARAMETER PROTOTYPE:        ALN */
 243         declare_name (repeat_pdb);
 244         get_cl_param(\
 245                             /*argc*/      argc          ,\
 246                             /*argv*/      argv          ,\
 247                             /*output*/    &le           ,\
 248                             /*Name*/      "-repeat_pdb"         ,\
 249                             /*Flag*/      &garbage      ,\
 250                             /*TYPE*/      "S"           ,\
 251                             /*OPTIONAL?*/ OPTIONAL      ,\
 252                             /*MAX Nval*/  1           ,\
 253                             /*DOC*/       "ND"          ,\
 254                             /*Parameter*/  &repeat_pdb,\
 255                             /*Def 1*/    "",\
 256                             /*Def 2*/     "stdin"       ,\
 257                             /*Min_value*/ "any"         ,\
 258                             /*Max Value*/ "any"          \
 259                    );
 260
 261 /*PARAMETER PROTOTYPE:    Nb to exclude     */
 262                get_cl_param(\
 263                             /*argc*/      argc          ,\
 264                             /*argv*/      argv          ,\
 265                             /*output*/    &le           ,\
 266                             /*Name*/      "-n_excluded_nb"    ,\
 267                             /*Flag*/      &garbage      ,\
 268                             /*TYPE*/      "D"           ,\
 269                             /*OPTIONAL?*/ OPTIONAL      ,\
 270                             /*MAX Nval*/  1             ,\
 271                             /*DOC*/       "Exclude the N Nb on each side of the central residue. -1 triggers an automatic setting equal to the window size corresponding to the sphere"          ,\
 272                             /*Parameter*/ &n_excluded_nb          ,\
 273                             /*Def 1*/    "-1"         ,\
 274                             /*Def 2*/    "1"             ,\
 275                             /*Min_value*/ "any"         ,\
 276                             /*Max Value*/ "any"          \
 277                    );
 278 /*PARAMETER PROTOTYPE:   diatances to count    */
 279                get_cl_param(\
 280                             /*argc*/      argc          ,\
 281                             /*argv*/      argv          ,\
 282                             /*output*/    &le           ,\
 283                             /*Name*/      "-similarity_threshold"    ,\
 284                             /*Flag*/      &garbage       ,\
 285                             /*TYPE*/      "F"            ,\
 286                             /*OPTIONAL?*/ OPTIONAL       ,\
 287                             /*MAX Nval*/  1              ,\
 288                             /*DOC*/       "ND"           ,\
 289                             /*Parameter*/ &similarity_threshold,\
 290                             /*Def 1*/    "70"             ,\
 291                             /*Def 2*/    "70"             ,\
 292                             /*Min_value*/ "any"          ,\
 293                             /*Max Value*/ "any"           \
 294                    );
 295 /*PARAMETER PROTOTYPE:   diatances to count    */
 296                get_cl_param(\
 297                             /*argc*/      argc          ,\
 298                             /*argv*/      argv          ,\
 299                             /*output*/    &le           ,\
 300                             /*Name*/      "-filter"    ,\
 301                             /*Flag*/      &garbage       ,\
 302                             /*TYPE*/      "F"            ,\
 303                             /*OPTIONAL?*/ OPTIONAL       ,\
 304                             /*MAX Nval*/  1              ,\
 305                             /*DOC*/       "Filter by only keeping the best quantile"           ,\
 306                             /*Parameter*/ &filter,\
 307                             /*Def 1*/    "1.00"             ,\
 308                             /*Def 2*/    "1.00"             ,\
 309                             /*Min_value*/ "-1.00"          ,\
 310                             /*Max Value*/ "1.00"           \
 311                    );
 312 /*PARAMETER PROTOTYPE:   diatances to count    */
 313                get_cl_param(\
 314                             /*argc*/      argc          ,\
 315                             /*argv*/      argv          ,\
 316                             /*output*/    &le           ,\
 317                             /*Name*/      "-filter_aln"    ,\
 318                             /*Flag*/      &garbage       ,\
 319                             /*TYPE*/      "D"            ,\
 320                             /*OPTIONAL?*/ OPTIONAL       ,\
 321                             /*MAX Nval*/  1              ,\
 322                             /*DOC*/       "Lower Case For Residues Filtered Out"           ,\
 323                             /*Parameter*/ &filter_aln,\
 324                             /*Def 1*/    "0"             ,\
 325                             /*Def 2*/    "1"             ,\
 326                             /*Min_value*/ "0"          ,\
 327                             /*Max Value*/ "1"           \
 328                    );
 329 /*PARAMETER PROTOTYPE:   diatances to count    */
 330                get_cl_param(\
 331                             /*argc*/      argc          ,\
 332                             /*argv*/      argv          ,\
 333                             /*output*/    &le           ,\
 334                             /*Name*/      "-irmsd_graph"    ,\
 335                             /*Flag*/      &garbage       ,\
 336                             /*TYPE*/      "D"            ,\
 337                             /*OPTIONAL?*/ OPTIONAL       ,\
 338                             /*MAX Nval*/  1              ,\
 339                             /*DOC*/       "Outputs the irmsd, position/position"           ,\
 340                             /*Parameter*/ &irmsd_graph,\
 341                             /*Def 1*/    "0"             ,\
 342                             /*Def 2*/    "1"             ,\
 343                             /*Min_value*/ "0"          ,\
 344                             /*Max Value*/ "1"           \
 345                    );
 346 /*PARAMETER PROTOTYPE:   diatances to count    */
 347                get_cl_param(\
 348                             /*argc*/      argc          ,\
 349                             /*argv*/      argv          ,\
 350                             /*output*/    &le           ,\
 351                             /*Name*/      "-nirmsd_graph"    ,\
 352                             /*Flag*/      &garbage       ,\
 353                             /*TYPE*/      "D"            ,\
 354                             /*OPTIONAL?*/ OPTIONAL       ,\
 355                             /*MAX Nval*/  1              ,\
 356                             /*DOC*/       "Outputs the NIRMSD VS N Removed Residues Curve"           ,\
 357                             /*Parameter*/ &nirmsd_graph,\
 358                             /*Def 1*/    "0"             ,\
 359                             /*Def 2*/    "1"             ,\
 360                             /*Min_value*/ "0"          ,\
 361                             /*Max Value*/ "1"           \
 362                    );
 363 /*PARAMETER PROTOTYPE:    -rmsd_threshold     */
 364                get_cl_param(\
 365                             /*argc*/      argc          ,\
 366                             /*argv*/      argv          ,\
 367                             /*output*/    &le           ,\
 368                             /*Name*/      "-md_threshold"     ,\
 369                             /*Flag*/      &garbage      ,\
 370                             /*TYPE*/      "F"           ,\
 371                             /*OPTIONAL?*/ OPTIONAL      ,\
 372                             /*MAX Nval*/  1             ,\
 373                             /*DOC*/       "ND"          ,\
 374                             /*Parameter*/ &md_threshold ,\
 375                             /*Def 1*/    "1"            ,\
 376                             /*Def 2*/    "1"             ,\
 377                             /*Min_value*/ "any"         ,\
 378                             /*Max Value*/ "any"         \
 379                    );
 380
 381 /*PARAMETER PROTOTYPE:    -maximum distances     */
 382                get_cl_param(\
 383                             /*argc*/      argc          ,\
 384                             /*argv*/      argv          ,\
 385                             /*output*/    &le           ,\
 386                             /*Name*/      "-maximum_distance"     ,\
 387                             /*Flag*/      &garbage      ,\
 388                             /*TYPE*/      "F"           ,\
 389                             /*OPTIONAL?*/ OPTIONAL      ,\
 390                             /*MAX Nval*/  1             ,\
 391                             /*DOC*/       "ND"          ,\
 392                             /*Parameter*/ &maximum_distance          ,\
 393                             /*Def 1*/    "10"            ,\
 394                             /*Def 2*/    "10"             ,\
 395                             /*Min_value*/ "any"         ,\
 396                             /*Max Value*/ "any"         \
 397                    );
 398
 399
 400 /*PARAMETER PROTOTYPE:    -print_rapdb     */
 401                get_cl_param(\
 402                             /*argc*/      argc          ,\
 403                             /*argv*/      argv          ,\
 404                             /*output*/    &le           ,\
 405                             /*Name*/      "-print_rapdb"     ,\
 406                             /*Flag*/      &garbage      ,\
 407                             /*TYPE*/      "D"           ,\
 408                             /*OPTIONAL?*/ OPTIONAL      ,\
 409                             /*MAX Nval*/  1             ,\
 410                             /*DOC*/       "Prints the neighborhood of each pair of aligned residues, along with the associated local score"          ,\
 411                             /*Parameter*/ &print_rapdb          ,\
 412                             /*Def 1*/    "0"            ,\
 413                             /*Def 2*/    "1"             ,\
 414                             /*Min_value*/ "any"         ,\
 415                             /*Max Value*/ "any"         \
 416                    );
 417
 418 /*PARAMETER PROTOTYPE:    RUN_NAME     */
 419                declare_name (run_name);
 420                get_cl_param(\
 421                             /*argc*/      argc          ,\
 422                             /*argv*/      argv          ,\
 423                             /*output*/    &le           ,\
 424                             /*Name*/      "-run_name"    ,\
 425                             /*Flag*/      &garbage      ,\
 426                             /*TYPE*/      "W_F"         ,\
 427                             /*OPTIONAL?*/ OPTIONAL      ,\
 428                             /*MAX Nval*/  1             ,\
 429                             /*DOC*/       "ND"          ,\
 430                             /*Parameter*/ &run_name      ,\
 431                             /*Def 1*/    "default"      ,\
 432                             /*Def 2*/    ""             ,\
 433                             /*Min_value*/ "default"         ,\
 434                             /*Max Value*/ "any"          \
 435                    );
 436 /*PARAMETER PROTOTYPE:    OUTFILE     */
 437 /*PARAMETER PROTOTYPE:    OUTFILE     */
 438                declare_name ( outfile);
 439                get_cl_param(\
 440                             /*argc*/      argc          ,\
 441                             /*argv*/      argv          ,\
 442                             /*output*/    &le           ,\
 443                             /*Name*/      "-outfile"    ,\
 444                             /*Flag*/      &garbage      ,\
 445                             /*TYPE*/      "W_F"         ,\
 446                             /*OPTIONAL?*/ OPTIONAL      ,\
 447                             /*MAX Nval*/  1             ,\
 448                             /*DOC*/       "ND"          ,\
 449                             /*Parameter*/ &outfile      ,\
 450                             /*Def 1*/    "no"      ,\
 451                             /*Def 2*/     "default"             ,\
 452                             /*Min_value*/ "default"         ,\
 453                             /*Max Value*/ "any"          \
 454                    );
 455 /*PARAMETER PROTOTYPE:    OUTFILE     */
 456                declare_name ( apdb_outfile);
 457                get_cl_param(\
 458                             /*argc*/      argc          ,\
 459                             /*argv*/      argv          ,\
 460                             /*output*/    &le           ,\
 461                             /*Name*/      "-apdb_outfile"    ,\
 462                             /*Flag*/      &garbage      ,\
 463                             /*TYPE*/      "W_F"         ,\
 464                             /*OPTIONAL?*/ OPTIONAL      ,\
 465                             /*MAX Nval*/  1             ,\
 466                             /*DOC*/       "ND"          ,\
 467                             /*Parameter*/ &apdb_outfile      ,\
 468                             /*Def 1*/    "stdout"      ,\
 469                             /*Def 2*/    "default"             ,\
 470                             /*Min_value*/ "any"         ,\
 471                             /*Max Value*/ "any"          \
 472                    );
 473
 474 /*PARAMETER PROTOTYPE:    OUTPUT_FORMAT    */
 475                out_aln_format=declare_char ( 200, STRING);
 476                n_out_aln_format=get_cl_param(\
 477                             /*argc*/      argc           ,\
 478                             /*argv*/      argv           ,\
 479                             /*output*/    &le            ,\
 480                             /*Name*/      "-output"      ,\
 481                             /*Flag*/      &garbage       ,\
 482                             /*TYPE*/      "S"            ,\
 483                             /*OPTIONAL?*/ OPTIONAL       ,\
 484                             /*MAX Nval*/  200              ,\
 485                             /*DOC*/       "ND"           ,\
 486                             /*Parameter*/ out_aln_format,\
 487                             /*Def 1*/    "score_html"           ,\
 488                             /*Def 2*/    ""             ,\
 489                             /*Min_value*/ "any"          ,\
 490                             /*Max Value*/ "any"           \
 491                    );
 492
 493
 494
 495 /*PARAMETER PROTOTYPE:    INFILE    */
 496                declare_name (color_mode);
 497                get_cl_param(\
 498                             /*argc*/      argc           ,\
 499                             /*argv*/      argv           ,\
 500                             /*output*/    &le            ,\
 501                             /*Name*/      "-color_mode"        ,\
 502                             /*Flag*/      &garbage       ,\
 503                             /*TYPE*/      "S"            ,\
 504                             /*OPTIONAL?*/ OPTIONAL       ,\
 505                             /*MAX Nval*/  1              ,\
 506                             /*DOC*/       "ND"           ,\
 507                             /*Parameter*/ &color_mode ,\
 508                             /*Def 1*/    "apdb"            ,\
 509                             /*Def 2*/    "irmsd"             ,\
 510                             /*Min_value*/ "any"           ,\
 511                             /*Max Value*/ "any"            \
 512                    );
 513 /*PARAMETER PROTOTYPE:    INFILE    */
 514                declare_name (output_res_num);
 515                get_cl_param(\
 516                             /*argc*/      argc           ,\
 517                             /*argv*/      argv           ,\
 518                             /*output*/    &le            ,\
 519                             /*Name*/      "-seqnos"        ,\
 520                             /*Flag*/      &garbage       ,\
 521                             /*TYPE*/      "S"            ,\
 522                             /*OPTIONAL?*/ OPTIONAL       ,\
 523                             /*MAX Nval*/  1              ,\
 524                             /*DOC*/       "ND"           ,\
 525                             /*Parameter*/ &output_res_num ,\
 526                             /*Def 1*/    "off"            ,\
 527                             /*Def 2*/    "on"             ,\
 528                             /*Min_value*/ "any"           ,\
 529                             /*Max Value*/ "any"            \
 530                    );
 531                declare_name (cache);
 532                get_cl_param(\
 533                             /*argc*/      argc          ,\
 534                             /*argv*/      argv          ,\
 535                             /*output*/    &le           ,\
 536                             /*Name*/      "-cache"    ,\
 537                             /*Flag*/      &garbage      ,\
 538                             /*TYPE*/      "W_F"         ,\
 539                             /*OPTIONAL?*/ OPTIONAL      ,\
 540                             /*MAX Nval*/  1             ,\
 541                             /*DOC*/       "use,ignore,update,local, directory name"          ,\
 542                             /*Parameter*/ &cache       ,\
 543                             /*Def 1*/    "use"      ,\
 544                             /*Def 2*/    "update"      ,\
 545                             /*Min_value*/ "any"         ,\
 546                             /*Max Value*/ "any"          \
 547                    );
 548
 549                 declare_name (local_mode);
 550                 get_cl_param(                                   \
 551                             /*argc*/      argc          ,\
 552                             /*argv*/      argv          ,\
 553                             /*output*/    &le           ,\
 554                             /*Name*/      "-local_mode"    ,\
 555                             /*Flag*/      &garbage      ,\
 556                             /*TYPE*/      "W_F"         ,\
 557                             /*OPTIONAL?*/ OPTIONAL      ,\
 558                             /*MAX Nval*/  1             ,\
 559                             /*DOC*/       "Mode for choosing the Neighborhood (bubble or window)\nWhen selecting window, maximum distance becomes the window 1/2 size, in residues\nWhen using sphere, maximum_distance is the sphere radius in Angstrom"          ,\
 560                             /*Parameter*/ &local_mode       ,\
 561                             /*Def 1*/    "sphere"      ,\
 562                             /*Def 2*/    "window"      ,\
 563                             /*Min_value*/ "any"         ,\
 564                             /*Max Value*/ "any"          \
 565                                           );
 566
 567 /*PARAMETER PROTOTYPE:        IN */
 568                 template_file_list=declare_char (100, STRING);
 569                 n_template_file=get_cl_param(                   \
 570                             /*argc*/      argc          , \
 571                             /*argv*/      argv          , \
 572                             /*output*/    &le           ,\
 573                             /*Name*/      "-template_file"         ,\
 574                             /*Flag*/      &garbage      ,\
 575                             /*TYPE*/      "S"           ,\
 576                             /*OPTIONAL?*/ OPTIONAL      ,\
 577                             /*MAX Nval*/  1000           ,\
 578                             /*DOC*/       "List of templates file for the sequences",\
 579                             /*Parameter*/ template_file_list     ,      \
 580                             /*Def 1*/    "_SELF_P_",\
 581                             /*Def 2*/     "stdin"       ,\
 582                             /*Min_value*/ "any"         ,\
 583                             /*Max Value*/ "any"          \
 584                                           );
 585                 /*PARAMETER PROTOTYPE:        MODE */
 586                 declare_name (mode);
 587                 get_cl_param(                             \
 588                             /*argc*/      argc          , \
 589                             /*argv*/      argv          , \
 590                             /*output*/    &le           ,\
 591                             /*Name*/      "-mode"         ,\
 592                             /*Flag*/      &garbage      ,\
 593                             /*TYPE*/      "S"           ,\
 594                             /*OPTIONAL?*/ OPTIONAL      ,\
 595                             /*MAX Nval*/  1           ,\
 596                             /*DOC*/       "Mode: irmsd, ",\
 597                             /*Parameter*/ &mode    ,    \
 598                             /*Def 1*/    "irmsd",\
 599                             /*Def 2*/     "stdin"       ,\
 600                             /*Min_value*/ "any"         ,\
 601                             /*Max Value*/ "any"          \
 602                                           );
 603
 604
 605
 606         get_cl_param(\
 607                             /*argc*/      argc             ,\
 608                             /*argv*/      argv             ,\
 609                             /*output*/    &le              ,\
 610                             /*Name*/      "-prot_min_sim"        ,\
 611                             /*Flag*/      &prot_min_sim        ,\
 612                             /*TYPE*/      "D"              ,\
 613                             /*OPTIONAL?*/ OPTIONAL         ,\
 614                             /*MAX Nval*/  1                ,\
 615                             /*DOC*/       "Minimum similarity between a sequence and its PDB target" ,\
 616                             /*Parameter*/ &prot_min_sim          ,\
 617                             /*Def 1*/     "0"             ,\
 618                             /*Def 2*/     "20"             ,\
 619                             /*Min_value*/ "any"            ,\
 620                             /*Max Value*/ "any"             \
 621                    );
 622  set_int_variable ("prot_min_sim", prot_min_sim);
 623
 624 get_cl_param(\
 625                             /*argc*/      argc             ,\
 626                             /*argv*/      argv             ,\
 627                             /*output*/    &le              ,\
 628                             /*Name*/      "-prot_max_sim"        ,\
 629                             /*Flag*/      &prot_max_sim        ,\
 630                             /*TYPE*/      "D"              ,\
 631                             /*OPTIONAL?*/ OPTIONAL         ,\
 632                             /*MAX Nval*/  1                ,\
 633                             /*DOC*/       "Maximum similarity between a sequence and its BLAST relatives" ,\
 634                             /*Parameter*/ &prot_max_sim          ,\
 635                             /*Def 1*/     "90"             ,\
 636                             /*Def 2*/     "100"             ,\
 637                             /*Min_value*/ "any"            ,\
 638                             /*Max Value*/ "any"             \
 639                    );
 640  set_int_variable ("prot_max_sim", prot_max_sim);
 641
 642 get_cl_param(\
 643                             /*argc*/      argc             ,\
 644                             /*argv*/      argv             ,\
 645                             /*output*/    &le              ,\
 646                             /*Name*/      "-prot_min_cov"        ,\
 647                             /*Flag*/      &prot_min_cov        ,\
 648                             /*TYPE*/      "D"              ,\
 649                             /*OPTIONAL?*/ OPTIONAL         ,\
 650                             /*MAX Nval*/  1                ,\
 651                             /*DOC*/       "Minimum coverage of a sequence by its BLAST relatives" ,\
 652                             /*Parameter*/ &prot_min_cov          ,\
 653                             /*Def 1*/     "0"             ,\
 654                             /*Def 2*/     "0"             ,\
 655                             /*Min_value*/ "any"            ,\
 656                             /*Max Value*/ "any"             \
 657                    );
 658 set_int_variable ("prot_min_cov", prot_min_cov);
 659
 660 get_cl_param(\
 661                             /*argc*/      argc             ,\
 662                             /*argv*/      argv             ,\
 663                             /*output*/    &le              ,\
 664                             /*Name*/      "-pdb_min_sim"        ,\
 665                             /*Flag*/      &pdb_min_sim        ,\
 666                             /*TYPE*/      "D"              ,\
 667                             /*OPTIONAL?*/ OPTIONAL         ,\
 668                             /*MAX Nval*/  1                ,\
 669                             /*DOC*/       "Minimum similarity between a sequence and its PDB target" ,\
 670                             /*Parameter*/ &pdb_min_sim          ,\
 671                             /*Def 1*/     "35"             ,\
 672                             /*Def 2*/     "35"             ,\
 673                             /*Min_value*/ "any"            ,\
 674                             /*Max Value*/ "any"             \
 675                    );
 676
 677  set_int_variable ("pdb_min_sim", pdb_min_sim);
 678  get_cl_param(                                                  \
 679                             /*argc*/      argc             ,\
 680                             /*argv*/      argv             ,\
 681                             /*output*/    &le              ,\
 682                             /*Name*/      "-pdb_max_sim"        ,\
 683                             /*Flag*/      &pdb_max_sim        ,\
 684                             /*TYPE*/      "D"              ,\
 685                             /*OPTIONAL?*/ OPTIONAL         ,\
 686                             /*MAX Nval*/  1                ,\
 687                             /*DOC*/       "Maximum similarity between a sequence and its PDB target" ,\
 688                             /*Parameter*/ &pdb_max_sim          ,\
 689                             /*Def 1*/     "100"             ,\
 690                             /*Def 2*/     "0"             ,\
 691                             /*Min_value*/ "any"            ,\
 692                             /*Max Value*/ "any"             \
 693                    );
 694  set_int_variable ("pdb_max_sim", pdb_max_sim);
 695  get_cl_param(                                                  \
 696                             /*argc*/      argc             ,\
 697                             /*argv*/      argv             ,\
 698                             /*output*/    &le              ,\
 699                             /*Name*/      "-pdb_min_cov"        ,\
 700                             /*Flag*/      &pdb_min_cov        ,\
 701                             /*TYPE*/      "D"              ,\
 702                             /*OPTIONAL?*/ OPTIONAL         ,\
 703                             /*MAX Nval*/  1                ,\
 704                             /*DOC*/       "Minimum coverage of a sequence by its PDB target" ,\
 705                             /*Parameter*/ &pdb_min_cov          ,\
 706                             /*Def 1*/     "50"             ,\
 707                             /*Def 2*/     "25"             ,\
 708                             /*Min_value*/ "any"            ,\
 709                             /*Max Value*/ "any"             \
 710                    );
 711 set_int_variable ("pdb_min_cov", pdb_min_cov);
 712
 713
 714
 715 declare_name (pdb_blast_server);
 716  get_cl_param(\
 717                             /*argc*/      argc          ,\
 718                             /*argv*/      argv          ,\
 719                             /*output*/    &le           ,\
 720                             /*Name*/      "-pdb_blast_server"    ,\
 721                             /*Flag*/      &garbage      ,\
 722                             /*TYPE*/      "W_F"         ,\
 723                             /*OPTIONAL?*/ OPTIONAL      ,\
 724                             /*MAX Nval*/  1             ,\
 725                             /*DOC*/       "ND"          ,\
 726                             /*Parameter*/&pdb_blast_server       ,\
 727                             /*Def 1*/    "EBI"      ,\
 728                             /*Def 2*/    "default"      ,\
 729                             /*Min_value*/ "any"         ,\
 730                             /*Max Value*/ "any"          \
 731                    );
 732 declare_name (prot_blast_server);
 733  get_cl_param(\
 734                             /*argc*/      argc          ,\
 735                             /*argv*/      argv          ,\
 736                             /*output*/    &le           ,\
 737                             /*Name*/      "-blast"    ,\
 738                             /*Flag*/      &garbage      ,\
 739                             /*TYPE*/      "W_F"         ,\
 740                             /*OPTIONAL?*/ OPTIONAL      ,\
 741                             /*MAX Nval*/  1             ,\
 742                             /*DOC*/       "ND"          ,\
 743                             /*Parameter*/&prot_blast_server       ,\
 744                             /*Def 1*/    ""      ,\
 745                             /*Def 2*/    ""      ,\
 746                             /*Min_value*/ "any"         ,\
 747                             /*Max Value*/ "any"          \
 748                    );
 749  //make sure that -blast and -blast_server are both supported blast>blast_server
 750  if ( !prot_blast_server[0])
 751    {
 752      get_cl_param(                                              \
 753                             /*argc*/      argc          ,\
 754                             /*argv*/      argv          ,\
 755                             /*output*/    &le           ,\
 756                             /*Name*/      "-blast_server"    ,\
 757                             /*Flag*/      &garbage      ,\
 758                             /*TYPE*/      "W_F"         ,\
 759                             /*OPTIONAL?*/ OPTIONAL      ,\
 760                             /*MAX Nval*/  1             ,\
 761                             /*DOC*/       "ND"          ,\
 762                             /*Parameter*/&prot_blast_server       ,\
 763                             /*Def 1*/    "EBI"      ,\
 764                             /*Def 2*/    "default"      ,\
 765                             /*Min_value*/ "any"         ,\
 766                             /*Max Value*/ "any"          \
 767                    );
 768    }
 769  set_string_variable ("blast_server", prot_blast_server);
 770
 771
 772
 773  declare_name (pdb_db);
 774  get_cl_param(\
 775                             /*argc*/      argc          ,\
 776                             /*argv*/      argv          ,\
 777                             /*output*/    &le           ,\
 778                             /*Name*/      "-pdb_db"    ,\
 779                             /*Flag*/      &garbage      ,\
 780                             /*TYPE*/      "W_F"         ,\
 781                             /*OPTIONAL?*/ OPTIONAL      ,\
 782                             /*MAX Nval*/  1             ,\
 783                             /*DOC*/       "Non Redundant PDB database"          ,\
 784                             /*Parameter*/&pdb_db       ,\
 785                             /*Def 1*/    "pdb"      ,\
 786                             /*Def 2*/    "default"      ,\
 787                             /*Min_value*/ "any"         ,\
 788                             /*Max Value*/ "any"          \
 789                    );
 790  set_string_variable ("pdb_db", pdb_db);
 791
 792
 793 declare_name (prot_db);
 794  get_cl_param(\
 795                             /*argc*/      argc          ,\
 796                             /*argv*/      argv          ,\
 797                             /*output*/    &le           ,\
 798                             /*Name*/      "-protein_db"    ,\
 799                             /*Flag*/      &garbage      ,\
 800                             /*TYPE*/      "W_F"         ,\
 801                             /*OPTIONAL?*/ OPTIONAL      ,\
 802                             /*MAX Nval*/  1             ,\
 803                             /*DOC*/       "ND"          ,\
 804                             /*Parameter*/&prot_db       ,\
 805                             /*Def 1*/    "uniprot"      ,\
 806                             /*Def 2*/    "default"      ,\
 807                             /*Min_value*/ "any"         ,\
 808                             /*Max Value*/ "any"          \
 809                    );
 810
 811  get_cl_param(                                                  \
 812                             /*argc*/      argc          ,\
 813                             /*argv*/      argv          ,\
 814                             /*output*/    &le           ,\
 815                             /*Name*/      "-gapped"    ,\
 816                             /*Flag*/      &garbage      ,\
 817                             /*TYPE*/      "D"         ,\
 818                             /*OPTIONAL?*/ OPTIONAL      ,\
 819                             /*MAX Nval*/  1             ,\
 820                             /*DOC*/       "ND"          ,\
 821                             /*Parameter*/&gapped       ,\
 822                             /*Def 1*/    "0"      ,\
 823                             /*Def 2*/    "1"      ,\
 824                             /*Min_value*/ "any"         ,\
 825                             /*Max Value*/ "any"          \
 826                    );
 827  get_cl_param(                                                  \
 828                             /*argc*/      argc          ,\
 829                             /*argv*/      argv          ,\
 830                             /*output*/    &le           ,\
 831                             /*Name*/      "-min_ncol"    ,\
 832                             /*Flag*/      &garbage      ,\
 833                             /*TYPE*/      "D"         ,\
 834                             /*OPTIONAL?*/ OPTIONAL      ,\
 835                             /*MAX Nval*/  1             ,\
 836                             /*DOC*/       "minimum number of columns (negative: fraction)"          ,\
 837                             /*Parameter*/&min_ncol       ,\
 838                             /*Def 1*/    "4"      ,\
 839                             /*Def 2*/    "1"      ,\
 840                             /*Min_value*/ "any"         ,\
 841                             /*Max Value*/ "any"          \
 842                    );
 843  // set the correct mode:
 844  if ( strm (argv[0], "trmsd"))sprintf (mode, "trmsd");
 845
 846  set_string_variable ("prot_db", prot_db);
 847
 848
 849                 if (argc==1){myexit (EXIT_SUCCESS);}
 850
 851                 if ( strm (outfile,"no"))n_out_aln_format=0;
 852
 853                 get_cl_param( argc, argv,&le, NULL,NULL,NULL,0,0,NULL);
 854                 prepare_cache (cache);
 855
 856
 857                 if (strm ( aln, ""))
 858                   sprintf ( aln, "%s", argv[1]);
 859
 860                 if (!is_aln (aln))
 861                   {
 862                     printf_exit (EXIT_FAILURE, stderr, "\n\n---- ERROR: File %s must be a valid alignment [FATAL:%s-%s]\n\n",aln,argv[0], PROGRAM);
 863                   }
 864
 865                 pdb_param=vcalloc ( 1, sizeof(Pdb_param));
 866
 867                 pdb_param->similarity_threshold=similarity_threshold;
 868
 869                 pdb_param->md_threshold=md_threshold;
 870                 pdb_param->maximum_distance=maximum_distance;
 871
 872                 if ( n_excluded_nb>0)
 873                   pdb_param->n_excluded_nb=n_excluded_nb;
 874                 else if ( n_excluded_nb==-1)
 875                   pdb_param->n_excluded_nb=(int)((float)maximum_distance/(float)1.57);
 876                 /* Exclude all the nb within the bubble at +1, +2, +n*/
 877                 pdb_param->print_rapdb=print_rapdb;
 878                 pdb_param->comparison_io=comparison_io;
 879
 880                 pdb_param->local_mode=local_mode;
 881                 pdb_param->color_mode=lower_string (color_mode);
 882                 pdb_param->filter=filter;
 883                 pdb_param->filter_aln=filter_aln;
 884                 pdb_param->irmsd_graph=irmsd_graph;
 885                 pdb_param->nirmsd_graph=nirmsd_graph;
 886
 887                 sprintf ( list_file[n_list++], "S%s", aln);
 888
 889
 890                 if (!strm (repeat_seq, ""))
 891                   {
 892
 893                     sprintf ( template_file_list[0], "%s", process_repeat (list_file[0], repeat_seq, repeat_pdb));
 894                     fprintf ( le, "\n##Turn a repeat List into a Template File\n");
 895                     le=display_file_content (le,template_file_list[0]);
 896                     fprintf ( le, "\n\n");
 897                   }
 898                 S=read_seq_in_n_list (list_file, n_list, NULL, NULL);
 899
 900                 le=display_sequences_names ( S,le,0, 0);
 901
 902                 if ( n_template_file)
 903                   {
 904                     fprintf ( le, "\nLooking For Sequence Templates:\n");
 905                     for ( a=0; a< n_template_file; a++)
 906                       {
 907                         fprintf ( le, "\n\tTemplate Type: [%s] Mode Or File: [%s] [Start", template_type2type_name(template_file_list[a]), template_file_list[a]);
 908                         S=seq2template_seq(S, template_file_list[a], F);
 909                         fprintf ( le, "]");
 910                       }
 911                   }
 912
 913                 if ( !strm (run_name, "default"))
 914                   {
 915                     F=parse_fname(run_name);
 916                     sprintf (F->name, "%s", F->full);
 917                   }
 918                 else
 919                   {
 920                     F=parse_fname (aln);
 921                   }
 922
 923                 for ( a=0; a< S->nseq; a++)
 924                   {
 925                     char *p;
 926
 927                     p=seq2T_value (S, a, "template_file", "_P_");
 928
 929                     if (p)sprintf (S->file[a], "%s",p);
 930                   }
 931
 932                 CL=declare_constraint_list ( S,NULL, NULL, 0,NULL, NULL);
 933                 CL->T=vcalloc (S->nseq,sizeof (Ca_trace*));
 934
 935
 936                 for ( n_pdb=0,a=0; a<S->nseq; a++)
 937                   {
 938                     if ( !is_pdb_file ( S->file[a])){CL->T[a]=NULL;continue;}
 939                     CL->T[a]=read_ca_trace    (S->file[a], "ATOM");
 940                     CL->T[a]=trim_ca_trace    (CL->T[a], S->seq[a]);
 941                     (CL->T[a])->pdb_param=pdb_param;
 942                     n_pdb++;
 943                   }
 944
 945                 A=declare_aln (S);
 946
 947
 948                 A->residue_case=KEEP_CASE;
 949                 A=main_read_aln(aln, A);
 950                 EA=copy_aln (A, EA);
 951                 A->CL=CL;
 952
 953                 if ( strm (apdb_outfile, "default"))
 954                   sprintf ( apdb_outfile, "%s.apdb_result", F->name);
 955
 956                 if ( n_pdb<2)
 957                   {
 958                     FILE *fp;
 959                     fp=vfopen (apdb_outfile, "w");
 960                     fprintf (fp, "\nYour Alignment Does Not Contain Enough Sequences With a known Structure\n");
 961                     fprintf (fp, "To Use APDB, your alignment must include at least TWO sequences with a known structure.\n");
 962                     fprintf (fp, "These sequences must be named according to their PDB identifier, followed by the chain index (if any) ex: 1fnkA\n");
 963                     fprintf (fp, "[FATAL:%s]\n", PROGRAM);
 964                     vfclose (fp);
 965                   }
 966                 else if ( strm (mode, "irmsd"))
 967                   {
 968                       EA=analyse_pdb ( A, EA, apdb_outfile);
 969                   }
 970                 else if ( strm (mode, "msa2tree") || strm (mode, "trmsd"))
 971                   {
 972                     EA=msa2struc_dist ( A, EA,F->name, gapped, min_ncol);
 973                   }
 974                 le=display_output_filename ( le, "APDB_RESULT", "APDB_RESULT_FORMAT_01", apdb_outfile, CHECK);
 975
 976                 if ( n_pdb>=2)
 977                   {
 978                     declare_name (file_name);
 979                     for ( a=0; a< n_out_aln_format; a++)
 980                       {
 981                         if ( strm2( outfile, "stdout", "stderr"))sprintf (file_name, "%s", outfile);
 982                         else if ( strm (outfile, "default"))
 983                           sprintf (file_name, "%s.%s",F->name, out_aln_format[a]);
 984                         else
 985                           sprintf (file_name, "%s.%s",outfile,out_aln_format[a]);
 986
 987                         output_format_aln (out_aln_format[a],A,EA,file_name);
 988                         le=display_output_filename ( le, "MSA", out_aln_format[a], file_name, CHECK);
 989                       }
 990                   }
 991                 return EXIT_SUCCESS;
 992     }
 993
 994
 995
 996 Constraint_list * set_constraint_list4align_pdb (Constraint_list *CL,int seq, char *dp_mode, char *local_mode, char *param_file)
 997 {
 998   static Constraint_list *PWCL;
 999   static Pdb_param *pdb_param;
1000   char **x;
1001   int n;
1002
1003   if ( !CL)
1004     {
1005       free_constraint_list (PWCL);
1006       return NULL;
1007     }
1008   else if ( !PWCL)
1009     {
1010       PWCL=declare_constraint_list ( CL->S,NULL, NULL, 0,NULL, NULL);
1011
1012       pdb_param=vcalloc ( 1, sizeof(Pdb_param));
1013       pdb_param->N_ca=0;
1014       pdb_param->max_delta=2.0;
1015       pdb_param->maximum_distance=14;
1016       declare_name (pdb_param->local_mode);
1017       sprintf (pdb_param->local_mode, "%s", local_mode);
1018       pdb_param->scale=50;
1019
1020       PWCL->pw_parameters_set=1;
1021       PWCL->S=CL->S;
1022       PWCL->lalign_n_top=10;
1023       PWCL->sw_min_dist=10;
1024
1025       PWCL->T=vcalloc ( (PWCL->S)->nseq, sizeof (Ca_trace*));
1026
1027       PWCL->extend_jit=0;
1028       PWCL->maximise=1;
1029       /*PWCL->gop=-40;*/
1030       PWCL->gop=-50;
1031       PWCL->gep=-20;
1032       sprintf (CL->matrix_for_aa_group, "vasiliky");
1033       PWCL->use_fragments=0;
1034       PWCL->ktup=0;
1035       PWCL->TG_MODE=1;
1036     }
1037
1038
1039    if ( param_file && check_file_exists ( param_file) )
1040         {
1041           if ( (x=get_parameter ( "-nca",              &n, param_file))!=NULL){pdb_param->N_ca=atoi(x[0]);free_char (x, -1);}
1042           if ( (x=get_parameter ( "-max_delta",        &n, param_file))!=NULL){pdb_param->max_delta=atof(x[0]);free_char (x, -1);}
1043           if ( (x=get_parameter ( "-maximum_distance", &n, param_file))!=NULL){pdb_param->maximum_distance=atoi(x[0]);free_char (x, -1);}
1044           if ( (x=get_parameter ( "-local_mode",       &n, param_file))!=NULL){sprintf (pdb_param->local_mode, "%s",x[0]);free_char (x, -1);}
1045           if ( (x=get_parameter ( "-scale",            &n, param_file))!=NULL){pdb_param->scale=atoi(x[0]);free_char (x, -1);}
1046           if ( (x=get_parameter ( "-gapopen", &n, param_file))!=NULL){PWCL->gop=atoi(x[0]);free_char (x, -1);}
1047           if ( (x=get_parameter ( "-gapext" , &n, param_file))!=NULL){PWCL->gep=atof(x[0]);free_char (x, -1);}
1048
1049         }
1050
1051
1052
1053
1054    sprintf ( PWCL->dp_mode, "%s", dp_mode);
1055
1056    if (strm (PWCL->dp_mode, "lalign"))sprintf (PWCL->dp_mode,"sim_pair_wise_lalign");
1057    else if (strm (PWCL->dp_mode, "sw"))sprintf (PWCL->dp_mode,"gotoh_pair_wise_sw");
1058
1059    local_mode=pdb_param->local_mode;
1060    if ( strm ( local_mode, "hasch_ca_trace_nb"))      PWCL->evaluate_residue_pair=evaluate_ca_trace_nb;
1061    else if ( strm ( local_mode, "hasch_ca_trace_bubble")) PWCL->evaluate_residue_pair=evaluate_ca_trace_bubble;
1062    else if ( strm ( local_mode, "hasch_ca_trace_sap1_bubble")) PWCL->evaluate_residue_pair=evaluate_ca_trace_sap1_bubble;
1063    else if ( strm ( local_mode, "hasch_ca_trace_sap2_bubble")) PWCL->evaluate_residue_pair=evaluate_ca_trace_sap2_bubble;
1064
1065    else if ( strm ( local_mode, "hasch_ca_trace_transversal")) PWCL->evaluate_residue_pair=evaluate_ca_trace_transversal;
1066    else if ( strm ( local_mode, "hasch_ca_trace_bubble_2")) PWCL->evaluate_residue_pair=evaluate_ca_trace_bubble_2;
1067    else if ( strm ( local_mode, "hasch_ca_trace_bubble_3")) PWCL->evaluate_residue_pair=evaluate_ca_trace_bubble_3;
1068    else if ( strm ( local_mode, "custom_pair_score_function1"))  PWCL->evaluate_residue_pair=custom_pair_score_function1;
1069    else if ( strm ( local_mode, "custom_pair_score_function2"))  PWCL->evaluate_residue_pair=custom_pair_score_function2;
1070    else if ( strm ( local_mode, "custom_pair_score_function3"))  PWCL->evaluate_residue_pair=custom_pair_score_function3;
1071    else if ( strm ( local_mode, "custom_pair_score_function4"))  PWCL->evaluate_residue_pair=custom_pair_score_function4;
1072    else if ( strm ( local_mode, "custom_pair_score_function5"))  PWCL->evaluate_residue_pair=custom_pair_score_function5;
1073    else if ( strm ( local_mode, "custom_pair_score_function6"))  PWCL->evaluate_residue_pair=custom_pair_score_function6;
1074    else if ( strm ( local_mode, "custom_pair_score_function7"))  PWCL->evaluate_residue_pair=custom_pair_score_function7;
1075    else if ( strm ( local_mode, "custom_pair_score_function8"))  PWCL->evaluate_residue_pair=custom_pair_score_function8;
1076    else if ( strm ( local_mode, "custom_pair_score_function9"))  PWCL->evaluate_residue_pair=custom_pair_score_function9;
1077    else if ( strm ( local_mode, "custom_pair_score_function10")) PWCL->evaluate_residue_pair=custom_pair_score_function10;
1078
1079
1080    else
1081      {
1082        fprintf ( stderr, "\n%s is an unknown hasch mode, [FATAL]\n", local_mode);
1083        myexit (EXIT_FAILURE);
1084      }
1085
1086    if ( PWCL->T[seq]);
1087    else
1088      {
1089        PWCL->T[seq]=read_ca_trace (is_pdb_struc((CL->S)->name[seq]), "ATOM");
1090        (PWCL->T[seq])->pdb_param=pdb_param;
1091        PWCL->T[seq]=trim_ca_trace (PWCL->T[seq], (CL->S)->seq[seq]);
1092        PWCL->T[seq]=hasch_ca_trace(PWCL->T[seq]);
1093
1094      }
1095
1096
1097   return PWCL;
1098 }
1099
1100
1101
1102 int evaluate_ca_trace_nb (Constraint_list *CL, int s1, int r1, int s2, int r2)
1103    {
1104
1105      return (int)(neighborhood_match(CL, s1,r1, s2, r2, (CL->T[s1])->Chain,(CL->T[s2])->Chain ));
1106    }
1107 int evaluate_ca_trace_sap2_bubble (Constraint_list *CL, int s1, int r1, int s2, int r2)
1108      {
1109
1110
1111
1112        return sap2_neighborhood_match (CL, s1, r1, s2, r2, (CL->T[s1])->Bubble,(CL->T[s2])->Bubble );
1113
1114      }
1115 int evaluate_ca_trace_sap1_bubble (Constraint_list *CL, int s1, int r1, int s2, int r2)
1116      {
1117        /*
1118          Function documentation: start
1119
1120          int evaluate_ca_trace_sap1_bubble (Constraint_list *CL, int s1, int s2, int r1, int r2)
1121          This function evaluates the cost for matching two residues:
1122
1123          a1 is the cost for matching the two neighborood ( bubble type), using sap
1124          a1: [0,+100], +100 is the best possible match.
1125          a2 is the residue type weight:
1126             min=worst substitution value
1127             best=best of r1/r1, r2/r2-min
1128
1129             a2=(r1/r2 -min)/best --> a1:[0, 100]
1130
1131          score=a1*a2-->[-inf, +10000];
1132        */
1133
1134
1135
1136        float a1;
1137
1138
1139        a1=(int) sap1_neighborhood_match (CL, s1, r1, s2, r2, (CL->T[s1])->Bubble,(CL->T[s2])->Bubble );
1140
1141        return (int)a1;
1142
1143
1144      }
1145 int evaluate_ca_trace_bubble (Constraint_list *CL, int s1, int r1, int s2, int r2)
1146      {
1147        /*
1148          Function documentation: start
1149
1150          int evaluate_ca_trace_bubble (Constraint_list *CL, int s1, int s2, int r1, int r2)
1151          This function evaluates the cost for matching two residues:
1152
1153          a1 is the cost for matching the two neighborood ( bubble type)
1154          a1: [-inf,+100-scale], +100-scale is the best possible match.
1155
1156        */
1157
1158
1159
1160        float a1;
1161
1162
1163
1164        a1=(int) neighborhood_match (CL, s1, r1, s2, r2, (CL->T[s1])->Bubble,(CL->T[s2])->Bubble )-((CL->T[s1])->pdb_param)->scale;
1165
1166               return a1;
1167
1168
1169      }
1170 int evaluate_ca_trace_transversal (Constraint_list *CL, int s1, int r1, int s2, int r2)
1171      {
1172        return (int)(transversal_match (CL, s1, r1, s2, r2, (CL->T[s1])->Transversal,(CL->T[s2])->Transversal ));
1173      }
1174
1175 int evaluate_ca_trace_bubble_3 (Constraint_list *CL, int s1, int r1, int s2, int r2)
1176      {
1177        /*This Mode evaluates :
1178
1179          1-The Bubble
1180          2-The Match of the transversal residues
1181        */
1182
1183        int a1, l1;
1184        int a2, l2;
1185        int a;
1186
1187        l1=MAX(( (CL->T[s1])->Chain )->nb[r1][0] ,((CL->T[s2])->Chain )->nb[r2][0]);
1188        l2=MAX(( (CL->T[s1])->Bubble)->nb[r1][0], ((CL->T[s2])->Bubble)->nb[r2][0]);
1189
1190        a1=(int)(neighborhood_match (CL, s1, r1, s2, r2, (CL->T[s1])->Bubble,(CL->T[s2])->Bubble ));
1191        a2=(int)(transversal_match (CL, s1, r1, s2, r2, (CL->T[s1])->Transversal,(CL->T[s2])->Transversal ));
1192
1193        if ( !l1 && !l2)return 0;
1194        a=(a1+a2)/2;
1195        return a;
1196      }
1197 int evaluate_ca_trace_bubble_2 (Constraint_list *CL, int s1, int r1, int s2, int r2)
1198      {
1199        /*This Mode evaluates :
1200          1-The Ca neighborhood
1201          2-The Bubble
1202        */
1203
1204
1205        return (int)((neighborhood_match (CL, s1, r1, s2, r2, (CL->T[s1])->Chain,(CL->T[s2])->Chain )));
1206      }
1207
1208
1209 /*********************************************************************************************/
1210 /*                                                                                           */
1211 /*         FUNCTIONS FOR COMPARING TWO NEIGHBORHOODS:START                                   */
1212 /*                                                                                           */
1213 /*********************************************************************************************/
1214 float matrix_match (Constraint_list *CL, int s1, int r1, int s2, int r2, Struct_nb *nbs1, Struct_nb *nbs2)
1215
1216      {
1217        /*
1218          Function documentation: start
1219
1220          float matrix_match (Constraint_list *CL, int s1, int s2, int r1, int r2, Struct_nb *nbs1, Struct_nb *nbs2)
1221          This function evaluates the matrix for matching two residues:
1222
1223             min=worst substitution value
1224             best=best of r1/r1, r2/r2-min
1225
1226             a2=(r1/r2 -min)/best --> a1:[0, 100]
1227
1228          score=a1*a2-->[-inf, +10000];
1229        */
1230
1231
1232
1233        float a2;
1234        float m1, m2, m;
1235        static float min=0;
1236        int a, b;
1237
1238        if ( !CL->M)
1239          {
1240            CL->M=read_matrice ( "pam250mt");
1241            min=CL->M[0][0];
1242            for ( a=0; a< 26; a++)
1243              for ( b=0; b< 26; b++)min=MIN(min, CL->M[a][b]);
1244          }
1245
1246        if ( r1<=0 || r2<=0)return 0;
1247        m1=CL->M[(CL->S)->seq[s1][r1-1]-'A'][(CL->S)->seq[s1][r1-1]-'A']-min;
1248        m2=CL->M[(CL->S)->seq[s2][r2-1]-'A'][(CL->S)->seq[s2][r2-1]-'A']-min;
1249        m=MAX(m1, m2);
1250        a2=(CL->M[(CL->S)->seq[s1][r1-1]-'A'][(CL->S)->seq[s2][r2-1]-'A']-min)/m;
1251
1252        return a2;
1253      }
1254
1255
1256 float transversal_match (Constraint_list *CL, int s1, int r1, int s2, int r2, Struct_nb *nbs1, Struct_nb *nbs2)
1257       {
1258         int a, l1, l2;
1259         float score=0;
1260         float delta, max_delta;
1261         float max;
1262         Pdb_param*PP;
1263
1264         PP=(CL->T[s1])->pdb_param;
1265         max_delta=PP->max_delta;
1266
1267         l1=nbs1->nb[r1][0];
1268         l2=nbs2->nb[r2][0];
1269
1270         if ( l1!=l2 || l1<(PP->N_ca)) return 0;
1271
1272
1273         max=MAX(l1, l2)*max_delta;
1274         for ( delta=0,a=0; a< l2 ; a++)
1275                {
1276
1277                    delta+=max_delta-FABS((nbs1->d_nb[r1][a]-nbs2->d_nb[r2][a]));
1278                }
1279         score=(delta*100)/max;
1280
1281
1282
1283        return score;
1284       }
1285
1286 float neighborhood_match (Constraint_list *CL, int s1, int r1, int s2, int r2, Struct_nb *nbs1, Struct_nb *nbs2)
1287       {
1288           static float **table;
1289           static int table_size;
1290           int a, b, l1, l2;
1291           float score=0;
1292           float ins, del, sub;
1293           float delta, max_delta;
1294           float max;
1295           Pdb_param*PP;
1296
1297
1298        PP=(CL->T[s1])->pdb_param;
1299        max_delta=PP->max_delta;
1300
1301
1302        if ( r1> 0 && r2 >0) {r1--; r2--;}
1303        else return 0;
1304
1305        l1=nbs1->nb[r1][0];
1306        l2=nbs2->nb[r2][0];
1307
1308        if (table_size< (MAX(l1, l2)+1))
1309          {
1310          table_size=MAX(l1, l2)+1;
1311          if ( table)free_float (table, -1);
1312          table=NULL;
1313          }
1314        if ( !table) table=declare_float (table_size, table_size);
1315
1316
1317        max=MAX(l1, l2)*max_delta;
1318        if ( max==0)return 0;
1319
1320
1321        table[0][0]=0;
1322        for ( b=1; b<=l2; b++)
1323            {
1324            table[0][b]=0;
1325            }
1326        for ( a=1; a<=l1; a++)
1327            {
1328            table[a][0]=0;
1329            for ( b=1; b<=l2 ; b++)
1330                {
1331
1332                    delta=max_delta-FABS((nbs1->d_nb[r1][a]-nbs2->d_nb[r2][b]));
1333
1334                    del=table[a-1][b];
1335                    ins=table[a][b-1];
1336                    sub= table[a-1][b-1]+delta;
1337
1338                    if ( del >= ins && del >= sub)score=del;
1339                    else if ( ins >= del && ins >= sub) score=ins;
1340                    else score=sub;
1341                    table[a][b]=score;
1342                }
1343            }
1344
1345
1346        score=((((score)*100)/max));
1347
1348
1349        return score;
1350       }
1351
1352 float sap1_neighborhood_match (Constraint_list *CL, int s1, int r1, int s2, int r2, Struct_nb *nbs1, Struct_nb *nbs2)
1353       {
1354         /*
1355           Function documentation: start
1356
1357           float sap1_neighborhood_match (Constraint_list *CL, int s1, int s2, int r1, int r2, Struct_nb *nbs1, Struct_nb *nbs2)
1358           This function is adapted from Taylor, Orengo, Protein Structure Alignment JMB 1989, (208)1-22
1359           It is the first function where
1360           score= A/(|dra-drb|+b)
1361
1362           Function documentation: end
1363         */
1364
1365           static float **table;
1366           static int table_size;
1367           int a, b, l1, l2;
1368           float score=0;
1369           float ins, del, sub;
1370           float delta;
1371           float max;
1372
1373           int A=50;
1374           int B=5;
1375
1376
1377
1378
1379
1380
1381        if ( r1> 0 && r2 >0) {r1--; r2--;}
1382        else return 0;
1383
1384        l1=nbs1->nb[r1][0];
1385        l2=nbs2->nb[r2][0];
1386
1387        if (table_size< (MAX(l1, l2)+1))
1388          {
1389          table_size=MAX(l1, l2)+1;
1390          if ( table)free_float (table, -1);
1391          table=NULL;
1392          }
1393        if ( !table) table=declare_float (table_size, table_size);
1394
1395
1396        max=MAX(l1, l2)*(A/B);
1397        if ( max==0)return 0;
1398
1399
1400        table[0][0]=0;
1401        for ( b=1; b<=l2; b++)
1402            {
1403            table[0][b]=0;
1404            }
1405        for ( a=1; a<=l1; a++)
1406            {
1407            table[a][0]=0;
1408            for ( b=1; b<=l2 ; b++)
1409                {
1410
1411                    delta=A/(FABS((nbs1->d_nb[r1][a]-nbs2->d_nb[r2][b]))+B);
1412
1413                    del=table[a-1][b];
1414                    ins=table[a][b-1];
1415                    sub= table[a-1][b-1]+delta;
1416                    if ( del >= ins && del >= sub)score=del;
1417                    else if ( ins >= del && ins >= sub) score=ins;
1418                    else score=sub;
1419                    table[a][b]=score;
1420                }
1421            }
1422
1423
1424        score=((score*100))/(max);
1425
1426
1427        return score;
1428       }
1429
1430 float sap2_neighborhood_match (Constraint_list *CL, int s1, int r1, int s2, int r2, Struct_nb *nbs1, Struct_nb *nbs2)
1431       {
1432         /*
1433           Function documentation: start
1434
1435           float sap1_neighborhood_match (Constraint_list *CL, int s1, int s2, int r1, int r2, Struct_nb *nbs1, Struct_nb *nbs2)
1436           This function is adapted from Taylor, Orengo, Protein Structure Alignment JMB 1989, (208)1-22
1437           It is the first function where
1438           score= A/(|dra-drb|+b)
1439
1440           Function documentation: end
1441         */
1442
1443           static float **table;
1444           static int table_size;
1445           int a, b, l1, l2;
1446           float score=0;
1447           float ins, del, sub;
1448           float delta;
1449           float max;
1450
1451           Amino_acid **pep1;
1452           Amino_acid **pep2;
1453           static Atom *vX_1, *vY_1, *vZ_1;
1454           static Atom *vX_2, *vY_2, *vZ_2;
1455           static Atom *ca1, *ca2;
1456           float val;
1457
1458           int A=50;
1459           int B=2;
1460
1461
1462
1463
1464        if ( r1> 0 && r2 >0) {r1--; r2--;}
1465        else return 0;
1466
1467        /*Make up the referencial*/
1468        pep1=(CL->T[s1])->peptide_chain;
1469        pep2=(CL->T[s2])->peptide_chain;
1470
1471        /*Get Referencial for CA1*/
1472        if ( (pep1[r1])->C)vX_1 =diff_atom(pep1[r1]->C,pep1[r1]->CA, vX_1);
1473        if ( (pep1[r1])->N)vY_1 =diff_atom(pep1[r1]->N,pep1[r1]->CA, vY_1);
1474        if ( (pep1[r1])->CB)vZ_1=diff_atom(pep1[r1]->CB,(pep1[r1])->CA,vZ_1);
1475        else vZ_1=add_atom (vX_1, vY_1, vZ_1);
1476
1477
1478
1479
1480
1481        /*Get Referencial for CA2*/
1482        if ( (pep2[r2])->C)vX_2 =diff_atom((pep2[r2])->C,(pep2[r2])->CA, vX_2);
1483        if ( (pep2[r2])->N)vY_2 =diff_atom((pep2[r2])->N,(pep2[r2])->CA, vY_2);
1484        if ( (pep2[r2])->CB)vZ_2=diff_atom((pep2[r2])->CB,(pep2[r2])->CA, vZ_2);
1485        else vZ_2=add_atom (vX_2, vY_2, vZ_2);
1486
1487
1488
1489
1490        /*END OF GETTING REFERENCIAL*/
1491
1492        /*Test
1493        if ( r1>1 && r2>1)
1494          {
1495          fprintf (stdout,"\n\t*******");
1496
1497          fprintf (stdout, "RESIDUE %d %c", r1, (CL->S)->seq[s1][r1]);
1498          if ( (pep1[r1])->CA)fprintf (stdout,"\n\tCA ");print_atom (pep1[r1]->CA );
1499          if ( (pep1[r1])->C)fprintf (stdout,"\n\tC  ");print_atom (pep1[r1]->C );
1500          if ( (pep1[r1])->N)fprintf (stdout,"\n\tN  ");print_atom (pep1[r1]->N );
1501          if ( (pep1[r1])->CB)fprintf (stdout,"\n\tCB ");print_atom (pep1[r1]->CB );
1502          fprintf (stdout,"\n\t*******");
1503          fprintf (stdout,"\n\tvX ");print_atom ( vX_1);
1504          fprintf (stdout,"\n\tvY ");print_atom ( vY_1);
1505          fprintf (stdout,"\n\tvZ ");print_atom ( vZ_1);
1506
1507          ca1= copy_atom ((pep1[r1-1])->CA, ca1);
1508          ca1 =diff_atom(ca1,(pep1[r1])->CA, ca1);
1509          fprintf (stdout,"\n\tca ");print_atom ( ca1);
1510          fprintf ( stdout, "\n\tSQ1=%d ", (int)square_atom(ca1));
1511          ca1=reframe_atom(vX_1, vY_1, vZ_1, ca1, ca1);
1512          fprintf ( stdout, "\n\tSQ2=%d ", (int)square_atom(ca1));
1513          fprintf (stdout,"\n\tca ");print_atom ( ca1);
1514          fprintf (stdout,"\n\n");
1515          }
1516        */
1517
1518        l1=nbs1->nb[r1][0];
1519        l2=nbs2->nb[r2][0];
1520
1521        if (table_size< (MAX(l1, l2)+1))
1522          {
1523          table_size=MAX(l1, l2)+1;
1524          if ( table)free_float (table, -1);
1525          table=NULL;
1526          }
1527        if ( !table) table=declare_float (table_size, table_size);
1528
1529
1530        max=MAX(l1, l2)*(A/B);
1531
1532        if ( max==0)return 0;
1533
1534
1535        table[0][0]=0;
1536        for ( b=1; b<=l2; b++)
1537            {
1538            table[0][b]=0;
1539            }
1540
1541        for ( a=1; a<=l1; a++)
1542            {
1543            ca1=copy_atom ((CL->T[s1])->structure[nbs1->nb[r1][a]], ca1);
1544            ca1=diff_atom(ca1,(pep1[r1])->CA, ca1);
1545            ca1=reframe_atom(vX_1, vY_1, vZ_1, ca1, ca1);
1546
1547            table[a][0]=0;
1548            for ( b=1; b<=l2 ; b++)
1549                {
1550                   ca2  =copy_atom((CL->T[s2])->structure[nbs2->nb[r2][b]], ca2);
1551                   ca2  =diff_atom(ca2,(pep2[r2])->CA, ca2);
1552                   ca2  =reframe_atom(vX_2, vY_2, vZ_2, ca2, ca2);
1553
1554                   ca2=diff_atom(ca2,ca1,ca2);
1555                   val=square_atom (ca2);
1556
1557                   val=(float)sqrt ((double)val);
1558
1559                   delta=A/(val+B);
1560
1561
1562                   del=table[a-1][b];
1563                   ins=table[a][b-1];
1564                   sub= table[a-1][b-1]+delta;
1565
1566                   if ( del >= ins && del >= sub)score=del;
1567                   else if ( ins >= del && ins >= sub) score=ins;
1568                   else score=sub;
1569                   table[a][b]=score;
1570                }
1571            }
1572
1573
1574        score=(((score*100))/(max)-50);
1575
1576
1577        return score;
1578       }
1579
1580 /*********************************************************************************************/
1581 /*                                                                                           */
1582 /*         APDB                                                                              */
1583 /*                                                                                           */
1584 /*********************************************************************************************/
1585 float **** irmsdmin_apdb_filtration ( Alignment *A, float ****residues, Constraint_list *CL, Pdb_param *PP, FILE *fp)
1586 {
1587   int s1, s2, a,col1, n,n2=0, t,flag;
1588   int **pos, **list;
1589   float nirmsd, min_nirmsd,max_nirmsd,ref_sum, sum, sum2;
1590   float **normalized_len;
1591
1592   normalized_len=declare_float (A->nseq+1, A->nseq+1);
1593   for (s1=0; s1<A->nseq; s1++)
1594     {
1595       int l1, l2, r1, r2, p;
1596       for (s2=0; s2<A->nseq; s2++)
1597         {
1598           for ( l1=l2=p=0; p< A->len_aln; p++)
1599             {
1600               r1=A->seq_al[s1][p];
1601               r2=A->seq_al[s2][p];
1602               if (!is_gap(r1) && isupper(r1))l1++;
1603               if (!is_gap(r2) && isupper(r2))l2++;
1604             }
1605           normalized_len[s1][s2]=MIN(l1,l2);
1606         }
1607     }
1608
1609   pos=aln2pos_simple (A, A->nseq);
1610   for ( s1=0; s1< A->nseq; s1++)
1611     for ( s2=0; s2<A->nseq; s2++)
1612       {
1613         if ( s1==s2) continue;
1614         else if (!(CL->T[A->order[s1][0]]) || !(CL->T[A->order[s2][0]]))continue;
1615
1616         list=declare_int (A->len_aln, 2);
1617
1618         for ( sum=0,n=0,col1=0; col1< A->len_aln; col1++)
1619           {
1620             if ( islower (A->seq_al[s1][col1]) || islower ( A->seq_al[s2][col1]))continue;
1621             else if ( pos[s1][col1]<=0 || pos[s2][col1]<=0 ) continue;
1622             else if ( residues[s1][s2][pos[s1][col1]-1][0]==0)continue;
1623
1624             list[n][0]=pos[s1][col1]-1;
1625             list[n][1]=(int)100000*residues[s1][s2][pos[s1][col1]-1][4];
1626             sum2+=residues[s1][s2][pos[s1][col1]-1][4];
1627             n++;
1628           }
1629
1630         if (n==0)return residues;
1631
1632         sort_int_inv (list, 2, 1,0, n-1);
1633         for (sum=0,a=0; a<n; a++)
1634           {
1635             sum+=list[a][1];
1636           }
1637         ref_sum=sum;
1638         nirmsd=min_nirmsd=max_nirmsd=sum/(n*n);
1639         t=0;
1640
1641
1642         /*1 Find the maximum*/
1643         sum=ref_sum;
1644         for (flag=0,a=0; a< n-1; a++)
1645           {
1646             sum-=list[a][1];
1647             nirmsd=sum/((n-(a+1))*(n-(a+1)));
1648             if (nirmsd<max_nirmsd)flag=1;
1649             if ((nirmsd>=max_nirmsd) && flag==1)break;
1650             n2=a;
1651           }
1652
1653         sum=ref_sum;
1654         for (a=0; a<n2-1; a++)
1655           {
1656             sum-=list[a][1];
1657             nirmsd=sum/((n-(a+1))*(n-(a+1)));
1658
1659
1660             if ( nirmsd<min_nirmsd)
1661               {
1662                 min_nirmsd=nirmsd;
1663                 t=a;
1664                 if ( PP->nirmsd_graph)
1665                   {
1666                     fprintf ( stdout, "\n_NIRMSD_GRAPH %s %s POS: %4d Removed: %4d NiRMSD: %.2f", A->name[s1], A->name[s2], list[a][0],a,(nirmsd/100000)*normalized_len[s1][s2]);
1667                   }
1668               }
1669           }
1670
1671         if ( PP->print_rapdb)
1672           {
1673             for ( a=0; a<n; a++)
1674               {
1675                 if      ( list[a][1]>0 && a<=t)fprintf ( stdout, "\nRAPDB QUANTILE REMOVE S1: %3d S2: %3d COL: %3d SCORE*100: %d", s1, s2, list[a][0], list[a][1]);
1676                 else if ( list[a][1]>0 && a>t)fprintf ( stdout, "\nRAPDB QUANTILE KEEP   S1: %3d S2: %3d COL: %3d SCORE*100: %d", s1, s2, list[a][0], list[a][1]);
1677               }
1678           }
1679
1680         fprintf ( stdout, "\n# MINIMISATION FILTER ON: NiRMSD minimsation resulted in the removal of %d [out of %d] Columns On the alignment %s Vs %s\n", t, n, A->name[s1], A->name[s2]);
1681         for ( a=0; a<=t; a++)
1682           {
1683
1684             residues[s1][s2][list[a][0]][0]=0;
1685             residues[s1][s2][list[a][0]][1]=0;
1686             residues[s1][s2][list[a][0]][2]=0;
1687             residues[s1][s2][list[a][0]][3]=0;
1688             residues[s1][s2][list[a][0]][4]=-1;
1689
1690           }
1691
1692         free_int (list, -1);
1693       }
1694   free_float (normalized_len, -1);
1695   return residues;
1696 }
1697 float **** quantile_apdb_filtration ( Alignment *A, float ****residues, Constraint_list *CL, Pdb_param *PP,FILE *fp)
1698 {
1699   int s1, s2, a,col1, n, t;
1700   int **pos, **list;
1701
1702   pos=aln2pos_simple (A, A->nseq);
1703   for ( s1=0; s1< A->nseq; s1++)
1704     for ( s2=0; s2<A->nseq; s2++)
1705       {
1706         if ( s1==s2) continue;
1707         else if (!(CL->T[A->order[s1][0]]) || !(CL->T[A->order[s2][0]]))continue;
1708
1709         list=declare_int (A->len_aln, 2);
1710
1711         for ( n=0,col1=0; col1< A->len_aln; col1++)
1712           {
1713             if ( islower (A->seq_al[s1][col1]) || islower ( A->seq_al[s2][col1]))continue;
1714             else if ( pos[s1][col1]<=0 || pos[s2][col1]<=0 ) continue;
1715
1716             list[n][0]=pos[s1][col1]-1;
1717             list[n][1]=(int)100*residues[s1][s2][pos[s1][col1]-1][4];
1718             n++;
1719
1720           }
1721
1722         sort_int_inv (list, 2, 1,0, n-1);
1723
1724         t=quantile_rank ( list,1, n,PP->filter);
1725
1726         if ( PP->print_rapdb)
1727           {
1728             for ( a=0; a<n; a++)
1729               {
1730                 if      ( list[a][1]>0 && a<t)fprintf ( stdout, "\nRAPDB QUANTILE REMOVE S1: %3d S2: %3d COL: %3d SCORE*100: %d", s1, s2, list[a][0], list[a][1]);
1731                 else if ( list[a][1]>0 && a>t)fprintf ( stdout, "\nRAPDB QUANTILE KEEP   S1: %3d S2: %3d COL: %3d SCORE*100: %d", s1, s2, list[a][0], list[a][1]);
1732               }
1733           }
1734
1735         for ( a=0; a<t; a++)
1736           {
1737
1738             residues[s1][s2][list[a][0]][0]=0;
1739             residues[s1][s2][list[a][0]][1]=0;
1740             residues[s1][s2][list[a][0]][2]=0;
1741             residues[s1][s2][list[a][0]][3]=0;
1742             residues[s1][s2][list[a][0]][4]=-1;
1743
1744           }
1745
1746         free_int (list, -1);
1747       }
1748
1749   return residues;
1750 }
1751 Alignment * analyse_pdb ( Alignment *A, Alignment *ST, char *results)
1752       {
1753         int s1,s2,r1, r2,b, p;
1754         int **pos;
1755         float **normalize_len;
1756         float m2, m4;
1757         float pair_tot=0, pair_m1, pair_m2, pair_m3, pair_m4, pair_m5, pair_len=0;
1758         float seq_tot, seq_m1, seq_m2, seq_m3, seq_m4, seq_m5,seq_len;
1759         float msa_tot, msa_m1, msa_m2, msa_m3, msa_m4, msa_m5, msa_len;
1760         float iRMSD_unit, iRMSD_max, iRMSD_min;
1761         float ****residues;
1762         Pdb_param *PP=NULL;
1763         Constraint_list *CL;
1764         char *average_file, *pairwise_file, *total_file, *irmsd_file=0;
1765         FILE *fp, *average,*pairwise, *total, *irmsd_graph=0;
1766
1767
1768         fp      =vfopen ( results, "w");
1769         pairwise=vfopen ((pairwise_file=vtmpnam (NULL)),"w");
1770         average =vfopen ((average_file =vtmpnam (NULL)),"w");
1771         total   =vfopen ((total_file   =vtmpnam (NULL)),"w");
1772
1773
1774         CL=A->CL;
1775
1776         for ( s1=0; s1< (A->S)->nseq; s1++)
1777           if ( CL->T[s1]){PP=(CL->T[s1])->pdb_param;break;}
1778
1779         if (PP->irmsd_graph)irmsd_graph   =vfopen ((irmsd_file =vtmpnam (NULL)),"w");
1780
1781         fprintf ( fp, "\nAPDB_RESULT_FORMAT_02\n");
1782         residues=analyse_pdb_residues ( A, A->CL,PP);
1783         if ( PP->filter>=0)residues=quantile_apdb_filtration (A, residues, A->CL,PP, fp);
1784         else if ( PP->filter<0)residues=irmsdmin_apdb_filtration (A, residues, A->CL,PP, fp);
1785
1786         pos=aln2pos_simple (A, A->nseq);
1787
1788
1789
1790
1791
1792         /*Compute the alignment length for normalization*/
1793           normalize_len=declare_float (A->nseq+1, A->nseq+1);
1794           for (s1=0; s1<A->nseq; s1++)
1795             {
1796               int l1, l2, r1, r2;
1797               for (s2=0; s2<A->nseq; s2++)
1798                 {
1799                   for ( l1=l2=p=0; p< A->len_aln; p++)
1800                     {
1801                       r1=A->seq_al[s1][p];
1802                       r2=A->seq_al[s2][p];
1803                       if (!is_gap(r1) && isupper(r1))l1++;
1804                       if (!is_gap(r2) && isupper(r2))l2++;
1805                     }
1806                   normalize_len[s1][s2]=MIN(l1,l2);
1807                 }
1808             }
1809
1810           msa_len=msa_tot=msa_m1=msa_m2=msa_m3=msa_m4=msa_m5=0;
1811
1812           for ( s1=0; s1< A->nseq; s1++)
1813             {
1814               if ( !(CL->T[A->order[s1][0]]))continue;
1815               seq_len=seq_tot=seq_m1=seq_m2=seq_m3=seq_m4=seq_m5=0;
1816               for ( s2=0; s2< A->nseq; s2++)
1817                 {
1818                   if ( s1==s2)continue;
1819                   if ( !(CL->T[A->order[s2][0]]))continue;
1820                   pair_tot=pair_m1=pair_m2=pair_m3=pair_m4=pair_m5=0;
1821                   for ( p=0; p< A->len_aln; p++)
1822                     {
1823                       r1=A->seq_al[s1][p];
1824                       r2=A->seq_al[s2][p];
1825                       b=pos[s1][p]-1;
1826
1827
1828                       if (PP->filter_aln)
1829                         {
1830                            if (is_gap(r1) || is_gap(r2) || residues[s1][s2][b][0]==0)
1831                              {
1832                                A->seq_al[s1][p]=tolower(r1);
1833                                A->seq_al[s2][p]=tolower(r2);
1834                              }
1835                            else
1836                              {
1837                                 A->seq_al[s1][p]=toupper(r1);
1838                                 A->seq_al[s2][p]=toupper(r2);
1839                              }
1840
1841                         }
1842
1843                       if ( PP->irmsd_graph && ( is_gap(r1) || is_gap(r2) || residues[s1][s2][b][0]==0))
1844                         {
1845
1846                           fprintf ( irmsd_graph, "\n_IRMSD_GRAPH %10s %10s ALN: %c%c iRMSD: -1.00", A->name[s1], A->name[s2],A->seq_al[s1][p], A->seq_al[s2][p]);
1847                         }
1848
1849                       if (is_gap(r1) || is_gap(r2) || residues[s1][s2][b][0]==0)continue;
1850                       pair_tot++;
1851
1852                       /*APDB*/
1853                       m2=(residues[s1][s2][b][2]*100)/residues[s1][s2][b][0];
1854                       if (m2>PP->similarity_threshold){pair_m3++;}
1855
1856                       /*iRMSD*/
1857
1858                       m4=residues[s1][s2][b][4];
1859
1860                       if ( PP->irmsd_graph )
1861                         {
1862                           fprintf ( irmsd_graph, "\nIRMSD_GRAPH %10s %10s ALN: %c%c iRMSD: %.2f", A->name[s1], A->name[s2],A->seq_al[s1][p], A->seq_al[s2][p], m4);
1863                           }
1864                       pair_m4+=m4;
1865                     }
1866                   pair_len=normalize_len[s1][s2];
1867                   if ( s1>s2)
1868                     {
1869
1870                       fprintf ( pairwise, "\n\n#PAIRWISE: %s Vs %s",A->name[s1], A->name[s2]);
1871                       fprintf ( pairwise, "\n\tPAIRWISE EVALUATED: %6.2f %%    [%s Vs %s] ",  (pair_len==0)?-1:(pair_tot*100)/pair_len,A->name[s1], A->name[s2]);
1872                       fprintf ( pairwise, "\n\tPAIRWISE APDB:      %6.2f %%    [%s Vs %s] ",  (pair_tot==0)?-1:(pair_m3*100)/pair_tot,A->name[s1], A->name[s2]);
1873                       fprintf ( pairwise, "\n\tPAIRWISE iRMSD:     %6.2f Angs [%s Vs %s]",  (pair_tot==0)?-1:pair_m4/pair_tot,A->name[s1], A->name[s2]);
1874                       fprintf ( pairwise, "\n\tPAIRWISE NiRMSD:    %6.2f Angs [%s Vs %s] [%d pos]", (pair_tot==0)?-1:(pair_m4*pair_len)/(pair_tot*pair_tot), A->name[s1], A->name[s2], (int)pair_tot);
1875                       fprintf ( pairwise, "\n\tRAPDB PAIRS PAIRWISE N_NONEMPTY_PAIRS %d N_MAXIMUM_PAIRS %d",(int) pair_tot, (int)pair_len);
1876                       msa_m3+=pair_m3;
1877                        msa_m4+=pair_m4;
1878                        msa_tot+=pair_tot;
1879                        msa_len+=pair_len;
1880                     }
1881                   seq_m3+=pair_m3;
1882                   seq_m4+=pair_m4;
1883                   seq_tot+=pair_tot;
1884                   seq_len+=pair_len;
1885
1886                 }
1887
1888               fprintf ( average, "\n\n#AVERAGE For Sequence %s", A->name[s1]);
1889               fprintf ( average, "\n\tAVERAGE  EVALUATED: %6.2f %%    [%s]", (seq_len==0)?-1:(seq_tot*100)/seq_len, A->name[s1]);
1890               fprintf ( average, "\n\tAVERAGE  APDB:      %6.2f %%    [%s]", (seq_tot==0)?-1:(seq_m3*100)/seq_tot, A->name[s1]);
1891               fprintf ( average, "\n\tAVERAGE  iRMSD:     %6.2f Angs [%s]", (seq_tot==0)?-1:seq_m4/seq_tot, A->name[s1]);
1892               fprintf ( average, "\n\tAVERAGE  NiRMS:     %6.2f Angs [%s]", (seq_tot==0)?-1:(seq_m4*seq_len)/(seq_tot*seq_tot), A->name[s1]);
1893               if ( strm (PP->color_mode, "apdb"))ST->score_seq[s1]=(seq_tot==0)?-1:(seq_m3*100)/pair_tot;
1894               if (PP->print_rapdb)fprintf (average, "\n\tRAPDB PAIRS AVERAGE N_NONEMPTY_PAIRS %d N_MAXIMUM_PAIRS %d", (int)pair_tot, (int)pair_len);
1895
1896               if ( strm (PP->color_mode, "irmsd"))ST->score_seq[s1]=(seq_tot==0)?-1:10*((seq_m4*pair_len)/(seq_tot*seq_tot));
1897
1898             }
1899           fprintf ( total, "\n\n#TOTAL for the Full MSA");
1900           fprintf ( total, "\n\tTOTAL     EVALUATED: %6.2f %%  ", (msa_len==0)?-1:(msa_tot*100)/msa_len);
1901           fprintf ( total, "\n\tTOTAL     APDB:      %6.2f %%  ", (msa_tot==0)?-1:(msa_m3*100)/msa_tot);
1902           fprintf ( total, "\n\tTOTAL     iRMSD:     %6.2f Angs", (msa_tot==0)?-1:msa_m4/msa_tot);
1903           fprintf ( total, "\n\tTOTAL     NiRMSD:    %6.2f Angs", (msa_tot==0)?-1:(msa_m4*msa_len)/(msa_tot*msa_tot));
1904           if (PP->print_rapdb)fprintf (total, "\n\tRAPDB PAIRS TOTAL N_NONEMPTY_PAIRS: %d N_MAXIMUM_PAIRS %d", (int)msa_tot, (int)msa_len);
1905
1906           if ( strm (PP->color_mode, "apdb")) ST->score_aln=ST->score=A->score_aln=A->score=(msa_tot==0)?-1:(msa_m3*100)/msa_tot;
1907           if ( strm (PP->color_mode, "irmsd"))ST->score_aln=ST->score=A->score_aln=A->score=(msa_tot==0)?-1:10*((msa_m4*msa_len)/(msa_tot*msa_tot));
1908
1909           vfclose (average);vfclose (total); vfclose (pairwise);if (PP->irmsd_graph)vfclose (irmsd_graph);
1910           fp=display_file_content (fp, pairwise_file);
1911           fp=display_file_content (fp, average_file);
1912           fp=display_file_content (fp, total_file);
1913           if ( PP->irmsd_graph)fp=display_file_content (fp, irmsd_file);
1914
1915           fprintf ( fp, "\n\n# EVALUATED: Fraction of Pairwise Columns Evaluated\n");
1916           fprintf ( fp, "# APDB:      Fraction of Correct Columns according to APDB\n");
1917           fprintf ( fp, "# iRMDS:     Average iRMSD over all evaluated columns\n");
1918           fprintf ( fp, "# NiRMDS:    iRMSD*MIN(L1,L2)/Number Evaluated Columns\n");
1919           fprintf ( fp, "# Main Parameter: -maximum_distance %.2f Angstrom\n", PP->maximum_distance);
1920
1921           fprintf ( fp, "# Undefined values are set to -1 and indicate LOW Alignment Quality\n");
1922           fp=print_program_information (fp, NULL);
1923
1924
1925
1926
1927           /*Color Output*/
1928           for (iRMSD_max=0,iRMSD_min=10000,s1=0; s1<A->nseq; s1++)
1929             for ( s2=0; s2< A->nseq; s2++)
1930               for (p=0; p<A->len_aln; p++)
1931                 {
1932                   if ( residues[s1][s2][p][4]>0)
1933                     {
1934                     iRMSD_max=MAX(iRMSD_max, residues[s1][s2][p][4]);
1935                     iRMSD_min=MAX(iRMSD_min, residues[s1][s2][p][4]);
1936                     }
1937
1938                 }
1939           iRMSD_unit=iRMSD_max/8;
1940
1941           for (p=0; p< A->len_aln; p++)
1942             for ( s1=0; s1< A->nseq; s1++)
1943               {
1944
1945                 for ( p=0; p< A->len_aln; p++)
1946                   {
1947                     r1=A->seq_al[s1][p];
1948                     b=pos[s1][p]-1;
1949                     if ( is_gap(r1) ||  !(CL->T[A->order[s1][0]]))
1950                       ST->seq_al[s1][p]=NO_COLOR_RESIDUE;
1951                     else
1952                       {
1953                         float tot_m2=0, tot_m4=0, v=0;
1954                         seq_m2=seq_m4=0;
1955
1956                         for (s2=0; s2< A->nseq; s2++)
1957                           {
1958                             r2=A->seq_al[s1][p];
1959                             if ( s1==s2) continue;
1960                             if (is_gap(r2) || !(CL->T[A->order[s1][0]]) || residues[s1][s2][b][0]==0)continue;
1961
1962                             seq_m2+=m2=(residues[s1][s2][b][2]*100)/residues[s1][s2][b][0];
1963                             tot_m2++;
1964
1965                             m4=residues[s1][s2][b][4];
1966                             if (m4>=0)
1967                               {
1968                                 seq_m4+=m4;
1969                                 tot_m4++;
1970                               }
1971                           }
1972
1973                         if (strm ( PP->color_mode, "apdb"))
1974                           {
1975                             if (tot_m2==0)v=NO_COLOR_RESIDUE;
1976                             else v=MIN((seq_m2/(10*tot_m2)),9);
1977                           }
1978                         else if ( strm (PP->color_mode, "irmsd"))
1979                           {
1980                             if ( tot_m4==0)v=NO_COLOR_RESIDUE;
1981                             else v=(8-(int)((seq_m4/(iRMSD_unit*tot_m4))))+1;
1982                           }
1983                         ST->seq_al[s1][p]=v;
1984
1985                       }
1986                   }
1987               }
1988           for ( p=0; p<A->len_aln; p++) ST->seq_al[A->nseq][p]=NO_COLOR_RESIDUE;
1989
1990
1991           ST->generic_comment=vcalloc ( 100, sizeof (int));
1992           if ( strm (PP->color_mode, "apdb"))
1993             {
1994               sprintf ( ST->generic_comment, "# APDB Evaluation: Color Range Blue-[0 %% -- 100 %%]-Red\n# Sequence Score: APDB\n# Local Score: APDB\n\n");
1995             }
1996           else if ( strm (PP->color_mode, "irmsd"))
1997             {
1998               sprintf ( ST->generic_comment, "\n# iRMSD Evaluation:\n# Sequence score: NiRMSD (Angstrom*10)\n# Local Score: iRMSD, Blue-[%.2f Ang. -- 0.00 Ang.]-Red \n", iRMSD_max);
1999             }
2000
2001           fprintf ( fp, "\n");
2002           vfclose (fp);
2003           free_int (pos, -1);
2004           return ST;
2005       }
2006 float **** analyse_pdb_residues ( Alignment *A, Constraint_list *CL, Pdb_param *pdb_param)
2007      {
2008
2009          int **pos;
2010          int s1, s2, rs1, rs2;
2011          int col1, col2;
2012          float ****distances;
2013
2014               /*Distances[Nseq][len_aln][4]
2015                 distances...[0]: Number of residues within the bubble
2016                 distances...[1]: Absolute difference of distance of residues within Bubble
2017                 distances...[2]: Number of residues within the bubble with Delta dist < md_threshold
2018                 distances ..[3]: Sum of squared difference of distances
2019                 distances ..[4]: iRMSD
2020               */
2021          float d1, d2,delta;
2022          int wd1, wd2;
2023          int in_bubble=0;
2024          int real_res1_col1=0;
2025          int real_res1_col2;
2026          int real_res2_col1;
2027          int real_res2_col2;
2028          Pdb_param *PP;
2029          int print_rapdb;
2030          float nrapdb, rapdb;
2031          Alignment *BA=NULL;
2032
2033          PP=pdb_param;
2034          print_rapdb=PP->print_rapdb;
2035
2036          distances=declare_arrayN(4, sizeof (float), A->nseq, A->nseq, 0, 0);
2037
2038          /*Pre-computation of the internal distances----> T[seq]->ca_dist[len][len]*/
2039          /*Can be avoided if distance_on_request set to 1 */
2040
2041          for ( s1=0; s1< A->nseq; s1++)
2042            {
2043              rs1=A->order[s1][0];
2044              if (CL->T[rs1] &&  !(CL->T[rs1])->ca_dist)(CL->T[rs1])->ca_dist=measure_ca_distances(CL->T[rs1]);
2045              for ( s2=0; s2< A->nseq; s2++)
2046                {
2047                  distances[s1][s2]=declare_float ( A->len_aln, 6);
2048                }
2049            }
2050          pos=aln2pos_simple (A, A->nseq);
2051
2052          for ( s1=0; s1< A->nseq; s1++)
2053            for ( col1=0; col1< A->len_aln; col1++)
2054              for ( s2=0; s2<A->nseq; s2++)
2055                {
2056                  rs1=A->order[s1][0];
2057                  rs2=A->order[s2][0];
2058                  rapdb=0;
2059                  nrapdb=0;
2060                  if ( s1==s2) continue;
2061                  else if (!(CL->T[rs1]) || !(CL->T[rs2]))continue;
2062                  else if ( islower (A->seq_al[s1][col1]) || islower ( A->seq_al[s2][col1]))continue;
2063                  else if ( pos[s1][col1]<=0 || pos[s2][col1]<=0 ) continue;
2064
2065                  if ( print_rapdb && s2>s1)
2066                    {
2067
2068                      fprintf ( stdout, "RAPDB S1: %s S2: %s POS %d %d %c %d %c ", A->name[s1], A->name[s2], col1+1, pos[s1][col1],A->seq_al[s1][col1], pos[s2][col1],A->seq_al[s2][col1]);
2069                      BA=copy_aln (A, BA);
2070                      lower_string (BA->seq_al[s1]);
2071                      lower_string (BA->seq_al[s2]);
2072                      BA->seq_al[s1][col1]=toupper (BA->seq_al[s1][col1]);
2073                      BA->seq_al[s2][col1]=toupper (BA->seq_al[s2][col1]);
2074                    }
2075
2076                  for ( col2=0; col2<A->len_aln; col2++)
2077                    {
2078
2079                      if (pos[s1][col2]<=0 || pos[s2][col2]<=0 )continue;
2080                      else if ( FABS((pos[s1][col2]-pos[s1][col1]))<=PP->n_excluded_nb)continue;
2081                      else if ( FABS((pos[s2][col2]-pos[s2][col1]))<=PP->n_excluded_nb)continue;
2082                      else if ( islower (A->seq_al[s1][col2]) || islower ( A->seq_al[s2][col2]))continue;
2083
2084                      real_res1_col1=pos[s1][col1]-1;
2085                      real_res1_col2=pos[s1][col2]-1;
2086
2087                      real_res2_col1=pos[s2][col1]-1;
2088                      real_res2_col2=pos[s2][col2]-1;
2089
2090                      d1=(CL->T[rs1])->ca_dist[real_res1_col1][real_res1_col2];
2091                      d2=(CL->T[rs2])->ca_dist[real_res2_col1][real_res2_col2];
2092
2093                      if ( d1==UNDEFINED || d2 == UNDEFINED) continue;
2094
2095
2096
2097                      if ( strm ( PP->local_mode, "sphere"))
2098                        {
2099                          in_bubble= (d1<PP->maximum_distance && d2<PP->maximum_distance)?1:0;              ;
2100                        }
2101                      else if ( strm ( PP->local_mode, "window"))
2102                        {
2103                          wd1=FABS((pos[s1][col2]-pos[s1][col1]));
2104                          wd2=FABS((pos[s2][col2]-pos[s2][col1]));
2105                          in_bubble= (wd1<PP->maximum_distance && wd2<PP->maximum_distance)?1:0;            ;
2106                        }
2107
2108                      if (in_bubble)
2109                        {
2110                          if ( print_rapdb && s2 >s1)
2111                            {
2112                              fprintf ( stdout, "NB %d %d %c %d %c ", col2, pos[s1][col2], A->seq_al[s1][col2], pos[s2][col2], A->seq_al[s2][col2]);
2113                              BA->seq_al[s1][col2]=toupper (BA->seq_al[s1][col2]);
2114                              BA->seq_al[s2][col2]=toupper (BA->seq_al[s2][col2]);
2115                            }
2116                          delta=FABS((d1-d2));
2117                          if (delta<PP->md_threshold)
2118                            distances[s1][s2][real_res1_col1][2]++;
2119                          distances[s1][s2][real_res1_col1][1]+=delta;
2120                          distances[s1][s2][real_res1_col1][0]++;
2121                          distances[s1][s2][real_res1_col1][3]+=delta*delta;
2122                          nrapdb++;
2123                          rapdb+=delta*delta;
2124                        }
2125                    }
2126
2127                  if ( nrapdb==0)distances[s1][s2][real_res1_col1][4]=-1;
2128                  else distances[s1][s2][real_res1_col1][4]=(float)sqrt((double)(rapdb/nrapdb));
2129
2130                  if ( print_rapdb && s2>s1)
2131                    {
2132                      if (nrapdb==0)
2133                        {
2134                          fprintf ( stdout, "APDB: UNDEFINED\n");
2135                        }
2136                      else
2137                        {
2138
2139                          fprintf ( stdout, " APDB: %.2f ",(float)sqrt((double)(rapdb/nrapdb)));
2140                          BA->residue_case=KEEP_CASE;unalign_residues (BA, s1, s2);
2141                          fprintf ( stdout,"SEQ1: %s %s SEQ2: %s %s\n", BA->name[s1], BA->seq_al[s1], BA->name[s2], BA->seq_al[s2]);
2142                        }
2143                    }
2144
2145                }
2146
2147          free_aln (BA);
2148          free_int (pos, -1);
2149          return distances;
2150      }
2151 int pair_res_suitable4trmsd    (int s1,int col1, int col2, Alignment *A, int **pos, Pdb_param *PP, Constraint_list *CL,int *s);
2152 int aln_column_contains_gap (Alignment *A, int c);
2153 float aln2ncol4trmsd(Alignment *A, int **pos, Constraint_list *CL, int **lc);
2154 int pair_columns_suitable4trmsd(int col1, int col2, Alignment *A, int **pos, Pdb_param *PP, Constraint_list *CL,int *s);
2155 int column_is_suitable4trmsd(int col1,Alignment *A, int **pos, Pdb_param *PP, Constraint_list *CL,int *s);
2156
2157
2158
2159 NT_node trmsdmat2tree (float **dm, int **count,Alignment *A);
2160 Alignment * msa2struc_dist ( Alignment *A, Alignment *ST, char *results, int gapped, int min_ncol4trmsd)
2161      {
2162
2163        int **pos, c;
2164        FILE *tl;
2165          int s1, s2, rs1, rs2;
2166          int col1, col2;
2167          float ****distances;
2168          float **dm,**tdm;
2169          int **count,**tcount;
2170          int print_subtrees=0;
2171          float min, max;
2172
2173               /*Distances[Nseq][len_aln][4]
2174                 distances...[0]: Number of residues within the bubble
2175                 distances...[1]: Absolute difference of distance of residues within Bubble
2176                 distances...[2]: Number of residues within the bubble with Delta dist < md_threshold
2177                 distances ..[3]: Sum of squared difference of distances
2178                 distances ..[4]: iRMSD
2179               */
2180          Pdb_param *pdb_param;
2181          Constraint_list *CL;
2182          int a, b, ncol, npos,n;
2183          float d1, d2,delta;
2184          int wd1, wd2;
2185          int in_bubble=0;
2186          int real_res1_col1=0;
2187          int real_res1_col2;
2188          int real_res2_col1;
2189          int real_res2_col2;
2190          Pdb_param *PP;
2191          int print_rapdb;
2192          float nrapdb, rapdb;
2193          Alignment *BA=NULL;
2194          NT_node *T0,*T1,*T2,*PT, *POS;
2195          NT_node BT0, BT10,BT50, BT100=NULL,RBT;
2196          char **pair_pos_list;
2197
2198          int ntree=0, ntree2;
2199
2200          Alignment *B;
2201          char *pos_list;
2202          char *tot_pos_list;
2203
2204          char *struc_tree10;
2205          char *struc_tree100;
2206          char *struc_tree50;
2207          char *struc_tree0;
2208          char *consense_file;
2209
2210          char *color_struc_tree;
2211          int **score;
2212          int proceed=1;
2213          int **lc;
2214          int used;
2215
2216          if (min_ncol4trmsd<0)
2217            {
2218              min_ncol4trmsd*=-1;
2219              min_ncol4trmsd=(min_ncol4trmsd*A->len_aln)/100;
2220            }
2221          else if ( min_ncol4trmsd>=A->len_aln)
2222            {
2223              min_ncol4trmsd=A->len_aln-1;
2224            }
2225
2226          lc=declare_int (A->nseq, 2);
2227          for (a=0; a<A->nseq; a++)lc[a][0]=a;
2228
2229          declare_name(tot_pos_list);
2230          sprintf ( tot_pos_list, "%s.struc_tree.list", results);
2231
2232          declare_name(consense_file);
2233          sprintf (consense_file, "%s.struc_tree.consense_output", results);
2234
2235          declare_name(pos_list);
2236          sprintf ( pos_list, "%s.pos_list", results);
2237
2238          declare_name(struc_tree0);
2239          sprintf ( struc_tree0, "%s.struc_tree.consensus",results);
2240
2241          declare_name(struc_tree10);
2242          sprintf ( struc_tree10, "%s.struc_tree10",results);
2243
2244          declare_name(struc_tree100);
2245          sprintf ( struc_tree100, "%s.struc_tree100",results);
2246
2247          declare_name(struc_tree50);
2248          sprintf ( struc_tree50, "%s.struc_tree50",results);
2249
2250          declare_name(color_struc_tree);
2251          sprintf ( color_struc_tree, "%s.struc_tree.html", results);
2252
2253          pair_pos_list=declare_char (A->len_aln*A->len_aln+1, 100);
2254          T1=vcalloc (A->len_aln*A->len_aln+1, sizeof (NT_node));
2255          T2=vcalloc (A->len_aln+1, sizeof (NT_node));
2256
2257          PT=vcalloc (A->len_aln*A->len_aln+1, sizeof (NT_node));
2258          POS=vcalloc (A->len_aln+1, sizeof (NT_node));
2259
2260          CL=A->CL;
2261
2262          //Check all sequences have a PDB structure
2263
2264          for (used=0,a=0; a<A->nseq; a++)
2265            {
2266              if ( ! seq2P_template_file(A->S,a))
2267                {
2268                  add_warning (stderr, "Sequence %s removed from the dataset [No Usable Structure]", A->name[a]);
2269                }
2270              else
2271                {
2272                  if (used!=a)
2273                    {
2274                      sprintf (A->name[used], "%s", A->name[a]);
2275                      sprintf (A->seq_al[used], "%s", A->seq_al[a]);
2276                      for (b=0; b<4; b++)A->order[used][b]=A->order[a][b];
2277                    }
2278                  used++;
2279                }
2280            }
2281
2282          A->nseq=used;
2283
2284          if (A->nseq<2)myexit (fprintf_error(stderr, "Two sequences at least must have a known structure"));
2285
2286          for ( s1=0; s1< (A->S)->nseq; s1++)
2287            if ( CL->T[s1]){PP=(CL->T[s1])->pdb_param;break;}
2288
2289          for ( s1=0; s1< A->nseq; s1++)
2290            {
2291              rs1=A->order[s1][0];
2292              if (CL->T[rs1] &&  !(CL->T[rs1])->ca_dist)(CL->T[rs1])->ca_dist=measure_ca_distances(CL->T[rs1]);
2293            }
2294          pos=aln2pos_simple (A, A->nseq);
2295
2296          dm=declare_float (A->nseq, A->nseq);
2297          count=declare_int (A->nseq, A->nseq);
2298          tdm=declare_float (A->nseq, A->nseq);
2299          tcount=declare_int (A->nseq, A->nseq);
2300
2301          PP->maximum_distance=1000;
2302          sprintf ( PP->local_mode, "sphere");
2303
2304          while ((npos=aln2ncol4trmsd(A,pos,CL,lc))<min_ncol4trmsd && A->nseq>1)
2305            {
2306
2307              sort_int_inv (lc,2, 1, 0,A->nseq-1);
2308              add_information (stderr, "Remove Sequence [%s] that contains %d un-suitable positions", A->name[lc[0][0]], lc[0][1]);
2309              A=remove_seq_from_aln (A, A->name[lc[0][0]]);
2310              ungap_aln (A);
2311              pos=aln2pos_simple (A, A->nseq);
2312            }
2313          if (!A->nseq)
2314            {
2315              myexit (fprintf_error(stderr,"No suitable pair of column supporting a tree"));
2316            }
2317          else
2318             fprintf ( stderr, "\n---- Number of usable positions: %d [%.2f %%]\n", npos, ((float)npos*100)/(float)A->len_aln);
2319
2320          tl=vfopen (tot_pos_list, "w");
2321          for (ncol=0,ntree=0, col1=0; col1< A->len_aln; col1++)
2322            {
2323              int w,tree, cont;
2324              //output_completion (stderr, col1, A->len_aln,1, "Sample Columns");
2325              if (!gapped && aln_column_contains_gap (A, col1))continue;
2326              for ( cont=1,ntree2=0,col2=0; col2<A->len_aln; col2++)
2327                {
2328                  for (s1=0; s1< A->nseq-1; s1++)
2329                    {
2330                      rs1=A->order[s1][0];
2331                      if (!pair_res_suitable4trmsd (s1,col1, col2, A, pos, PP, CL, &w))continue;
2332                      for ( s2=s1+1; s2<A->nseq; s2++)
2333                        {
2334                          if (!pair_res_suitable4trmsd (s2,col1, col2, A, pos, PP, CL, &w))continue;
2335
2336                          rs2=A->order[s2][0];
2337                          real_res1_col1=pos[s1][col1]-1;
2338                          real_res1_col2=pos[s1][col2]-1;
2339                          real_res2_col1=pos[s2][col1]-1;
2340                          real_res2_col2=pos[s2][col2]-1;
2341
2342                          d1=(CL->T[rs1])->ca_dist[real_res1_col1][real_res1_col2];
2343                          d2=(CL->T[rs2])->ca_dist[real_res2_col1][real_res2_col2];
2344
2345                          delta=FABS((d1-d2));
2346                          dm[s1][s2]=dm[s2][s1]+=delta;
2347                          tdm[s1][s2]=tdm[s2][s1]+=delta;
2348                          tcount[s1][s2]++;
2349                          tcount[s2][s1]++;
2350
2351                          count[s1][s2]++;
2352                          count[s2][s1]++;
2353                        }
2354                    }
2355                }
2356
2357
2358
2359              if ((POS[col1]=trmsdmat2tree (dm, count, A)))
2360                {
2361                  T1[ntree]=POS[col1];
2362                  fprintf (tl, "\n>Tree_%d Column\n", col1+1);
2363                  print_tree (T1[ntree], "newick", tl);
2364                  ntree++;
2365                }
2366            }
2367          vfclose (tl);
2368
2369          if (!ntree){fprintf ( stderr, "\nERROR: No suitable pair of column supporting a tree [FATAL]\n"); exit (EXIT_SUCCESS);}
2370
2371          score=treelist2avg_treecmp (T1, NULL);
2372          display_output_filename( stderr,"TreeList","newick",tot_pos_list, CHECK);
2373
2374          if (treelist_file2consense (tot_pos_list, NULL, consense_file))
2375            {
2376              display_output_filename( stderr,"ConsenseTree","phylip",consense_file, CHECK);
2377            }
2378          else
2379            {
2380              fprintf ( stderr, "\nPhylip is not installed: the program could not produce the consense output. This is not mandatory but useful");
2381            }
2382
2383          //consensus tree
2384
2385          if ((BT100=treelist2filtered_bootstrap (T1, NULL,score, 1.0)))
2386                {
2387                  vfclose (print_tree (BT100,"newick", vfopen (struc_tree0, "w")));
2388                  display_output_filename( stderr,"Tree","newick",struc_tree0, CHECK);
2389                }
2390          if (print_subtrees)
2391            {
2392
2393              if ( (BT0=trmsdmat2tree (tdm, tcount, A)))
2394                {
2395                  vfclose (print_tree (BT0,"newick", vfopen (struc_tree0, "w")));
2396                  display_output_filename( stderr,"Tree","newick",struc_tree0, CHECK);
2397                }
2398              if ((BT10=treelist2filtered_bootstrap (T1, NULL,score, 0.1)))
2399                {
2400                  vfclose (print_tree (BT10,"newick", vfopen (struc_tree10, "w")));
2401                  display_output_filename( stderr,"Tree","newick",struc_tree10, CHECK);
2402                }
2403
2404              if ((BT50=treelist2filtered_bootstrap (T1, NULL, score,0.5)))
2405                {
2406                  vfclose (print_tree (BT50,"newick", vfopen (struc_tree50, "w")));
2407                  display_output_filename( stderr,"Tree","newick",struc_tree50, CHECK);
2408                }
2409            }
2410
2411
2412          if (!BT100)BT100=treelist2filtered_bootstrap (T1, NULL,score, 1.0);
2413
2414          RBT=BT100;
2415          if (RBT)
2416            {
2417              B=copy_aln (A, NULL);
2418              for (a=0; a<A->len_aln; a++)
2419                {
2420                  int score;
2421                  Tree_sim *S=NULL;
2422
2423                  if (POS[a])
2424                    {
2425                      S=tree_cmp (POS[a], RBT);
2426                      score=S->uw/10;
2427                    }
2428                  else
2429                    {
2430                      score=NO_COLOR_RESIDUE;
2431                    }
2432
2433                  for (b=0; b<B->nseq; b++)
2434                    {
2435                      if ( is_gap (B->seq_al[b][a]) || score == NO_COLOR_RESIDUE)
2436                        {
2437                          B->seq_al[b][a]=NO_COLOR_RESIDUE;
2438                        }
2439                      else
2440                        {
2441                          B->seq_al[b][a]=S->uw/10;
2442                        }
2443                    }
2444                  if (S)vfree (S);
2445                }
2446
2447              output_format_aln ("score_html", A,B,color_struc_tree);
2448              display_output_filename( stderr,"Colored MSA","score_html",color_struc_tree, CHECK);
2449              free_aln (BA);
2450              fprintf ( stderr, "\n");
2451            }
2452          fprintf ( stderr, "\n");
2453          free_int (pos, -1);
2454          exit (EXIT_SUCCESS);
2455          return NULL;
2456      }
2457 NT_node trmsdmat2tree (float **dm, int **count,Alignment *A)
2458 {
2459   float min, max;
2460   int s1, s2;
2461   NT_node T;
2462   int ns;
2463   int **dm_int;
2464
2465   ns=A->nseq;
2466   for (s1=0; s1<ns-1; s1++)
2467     for (s2=s1+1; s2<ns; s2++)
2468       {
2469         if ( count [s1][s2])dm[s1][s2]=dm[s2][s1]=dm[s1][s2]/(float)count[s1][s2];
2470         else
2471           {
2472             return NULL;
2473           }
2474         if (s1==0 && s2==1)min=max=dm[s1][s2];
2475         min=MIN(dm[s1][s2], min);
2476         max=MAX(dm[s1][s2], max);
2477       }
2478   dm_int=declare_int (ns, ns);
2479   for (s1=0; s1<A->nseq-1; s1++)
2480     for (s2=s1+1; s2<A->nseq; s2++)
2481       {
2482         dm_int[s1][s2]=dm_int[s2][s1]=((dm[s1][s2])/(max))*100;
2483       }
2484   T=compute_std_tree_2(A, dm_int, "_TMODE_upgma");
2485   free_int (dm_int, -1);
2486   for (s1=0; s1<ns; s1++)for ( s2=0; s2<ns; s2++){dm[s1][s2]=count[s1][s2]=0;}
2487   return T;
2488 }
2489
2490 int pair_res_suitable4trmsd    (int s1,int col1, int col2, Alignment *A, int **pos, Pdb_param *PP, Constraint_list *CL,int *s)
2491 {
2492    int rs;
2493    rs=A->order[s1][0];
2494    if ( !(CL->T[rs])){s[0]=s1; return 0;}
2495    else if (is_gap (A->seq_al[s1][col1])){s[0]=s1;return 0;}
2496    else if (is_gap (A->seq_al[s1][col2])){s[0]=s1;return 0;}
2497
2498    else if (islower(A->seq_al[s1][col1])){s[0]=s1; return 0;}
2499    else if (islower(A->seq_al[s1][col2])){s[0]=s1; return 0;}
2500
2501    else if ( FABS(((pos[s1][col2])-(pos[s1][col1])))<=PP->n_excluded_nb){s[0]=s1;return 0;}
2502    else if ((CL->T[rs])->ca_dist[pos[s1][col1]-1][pos[s1][col2]-1]==UNDEFINED){s[0]=s1;return 0;}
2503    return 1;
2504 }
2505 int pair_columns_suitable4trmsd(int col1, int col2, Alignment *A, int **pos, Pdb_param *PP, Constraint_list *CL,int *s)
2506 {
2507   int s1;
2508   if (!column_is_suitable4trmsd (col1, A, pos, PP, CL,s))return 0;
2509   if (!column_is_suitable4trmsd (col2, A, pos, PP, CL,s))return 0;
2510   for (s1=0; s1<A->nseq; s1++)
2511     {
2512       int rs, rr1, rr2;
2513
2514       rs=A->order[s1][0];
2515       if ( FABS(((pos[s1][col2])-(pos[s1][col1])))<=PP->n_excluded_nb){s[0]=s1;return 0;}
2516       if ((CL->T[rs])->ca_dist[pos[s1][col1]-1][pos[s1][col2]-1]==UNDEFINED){s[0]=s1;return 0;}
2517       rr1=pos[s1][col1]-1;
2518       rr2=pos[s1][col2]-1;
2519       if ((CL->T[rs])->ca_dist[rr1][rr2]>PP->maximum_distance){s[0]=s1;return 0;}
2520     }
2521   return 1;
2522 }
2523 int column_is_suitable4trmsd(int col1,Alignment *A, int **pos, Pdb_param *PP, Constraint_list *CL,int *s)
2524 {
2525   int s1;
2526   for ( s1=0; s1<A->nseq; s1++)
2527     {
2528       int rs;
2529       rs=A->order[s1][0];
2530       if ( !(CL->T[rs])){s[0]=s1; return 0;}
2531       else if (is_gap (A->seq_al[s1][col1])){s[0]=s1;return 0;}
2532       else if (islower(A->seq_al[s1][col1])){s[0]=s1; return 0;}
2533     }
2534   return 1;
2535 }
2536 int aln_column_contains_gap (Alignment *A, int c)
2537 {
2538   int a, b;
2539   if ( !A || c>=A->len_aln || c<0)
2540     {
2541       printf ( "\nERROR: values out of range in aln_column_contains_gap [FATL:%s]\n", PROGRAM);
2542       exit (EXIT_FAILURE);
2543     }
2544   for ( a=0; a<A->nseq; a++) if ( is_gap(A->seq_al[a][c]))return 1;
2545   return 0;
2546 }
2547
2548
2549 float aln2ncol4trmsd(Alignment *A, int **pos, Constraint_list *CL, int **lc)
2550 {
2551   //This function estimates the number of columns suitable for constructing a trmsd
2552   int col1, s1, ncol, n, rs1, real_res1_col1;
2553
2554   for (s1=0; s1<A->nseq; s1++){lc[s1][0]=s1; lc[s1][1]=0;}
2555   for (ncol=0,col1=0; col1< A->len_aln; col1++)
2556     {
2557       for (n=0,s1=0; s1<A->nseq; s1++)
2558         {
2559           real_res1_col1=pos[s1][col1]-1;
2560           rs1=A->order[s1][0];
2561
2562           if (real_res1_col1<0)lc[s1][1]++;
2563           else if (!((CL->T[A->order[s1][0]])->ca[real_res1_col1]))lc[s1][1]++;
2564           else n++;
2565         }
2566       if (n==A->nseq)
2567         {
2568           ncol++;
2569         }
2570     }
2571   return ncol;
2572 }
2573
2574 float square_atom ( Atom *X)
2575 {
2576
2577   return X->x*X->x + X->y*X->y + X->z*X->z;
2578 }
2579 Atom* reframe_atom ( Atom *X, Atom*Y, Atom *Z, Atom *IN, Atom *R)
2580      {
2581        float new_x, new_y, new_z;
2582
2583        if ( R==NULL)R=vcalloc ( 1, sizeof (Atom));
2584
2585
2586         new_x= X->x*IN->x + Y->x*IN->y +Z->x*IN->z;
2587         new_y= X->y*IN->x + Y->y*IN->y +Z->y*IN->z;
2588         new_z= X->z*IN->x + Y->z*IN->y +Z->z*IN->z;
2589
2590         R->x=new_x;
2591         R->y=new_y;
2592         R->z=new_z;
2593        return R;
2594      }
2595
2596 Atom* add_atom ( Atom *A, Atom*B, Atom *R)
2597 {
2598   if ( R==NULL)R=vcalloc ( 1, sizeof (Atom));
2599
2600   R->x=A->x+B->x;
2601   R->y=A->y+B->y;
2602   R->z=A->z+B->z;
2603
2604   return R;
2605 }
2606 Atom* diff_atom ( Atom *A, Atom*B, Atom *R)
2607 {
2608   if ( R==NULL)R=vcalloc ( 1, sizeof (Atom));
2609
2610   R->x=A->x-B->x;
2611   R->y=A->y-B->y;
2612   R->z=A->z-B->z;
2613
2614   return R;
2615 }
2616
2617 Atom * copy_atom ( Atom *A, Atom*R)
2618 {
2619   if ( R==NULL)R=vcalloc ( 1, sizeof (Atom));
2620   R->num=A->num;
2621   R->res_num=A->res_num;
2622   R->x=A->x;
2623   R->y=A->y;
2624   R->z=A->z;
2625
2626   sprintf( R->type, "%s", A->type);
2627   return R;
2628 }
2629  void print_atom (Atom *A)
2630 {
2631   fprintf ( stdout, "%.2f %.2f %.2f", A->x, A->y, A->z);
2632 }
2633 /************************************************************************/
2634 /*                                                                      */
2635 /*                            NUSSINOV                                  */
2636 /*                                                                      */
2637 /************************************************************************/
2638
2639 /*---------prototypes ----------*/
2640 static void computeBasePairMatrix(int**M,char*S,int l, int T);
2641 static int backtrack(int a,int b,int**M,char*S,char*P, int T);
2642
2643
2644
2645 static int basePair(char x, char y)
2646 {
2647   static short **mat;
2648
2649   if (!mat)
2650     {
2651       char alp[20];
2652       int a, b, c1, c2, lc1, lc2;
2653       mat=declare_short (256, 256);
2654       sprintf ( alp, "AGCTUagctu");
2655       for (a=0; a<strlen (alp); a++)
2656         {
2657           for (b=a; b<strlen (alp)-1; b++)
2658             {
2659               c1=alp[a];c2=alp[b];
2660               lc1=tolower(c1); lc2=tolower(c2);
2661               if ( lc1=='g' && lc2=='c')
2662                 mat[c1][c2]=1;
2663               else if ( lc1=='a' && lc2=='u')
2664                 mat[c1][c2]=1;
2665               else if ( lc1=='u' && lc2=='g')
2666                 mat [c1][c2]=1;
2667               mat[c2][c1]=mat[c1][c2];
2668             }
2669         }
2670     }
2671   return (int)mat[(int)x][(int)y];
2672 }
2673
2674
2675
2676 /* ------------------------------------------------------------ */
2677
2678 char *nussinov(char *S, int THRESHOLD)
2679 {
2680   char *paren;
2681   int i;
2682
2683   /*-------------------------------
2684     S is RNA sequence
2685     paren is parenthesis expression for
2686     optimal RNA secondary structure
2687     THRESHOLD: Min distance between two paired residues
2688     -------------------------------*/
2689
2690   int **numBasePairs;
2691    int n;
2692
2693    /*----- initialization  --*/
2694    n = strlen(S);
2695    paren=vcalloc (n+1, sizeof (char));
2696    numBasePairs=declare_int (n,n);
2697
2698    for (i=0;i<n;i++) paren[i]='.';
2699    paren[n]='\0'; // paren is string of same length as S
2700    computeBasePairMatrix(numBasePairs,S,n, THRESHOLD);
2701    backtrack(0,n-1,numBasePairs,S,paren, THRESHOLD);
2702    free_int (numBasePairs, -1);
2703    return paren;
2704 }
2705
2706 static void computeBasePairMatrix(int** numBasePairs,char *S,int n, int THRESHOLD)
2707 {
2708    int i,j,d,k,max,val,index;
2709
2710    for (d = THRESHOLD; d < n; d++){
2711      for(i=0; i < n; i++)
2712        {
2713         j=i+d;
2714         if (j < n){
2715           max=0;
2716           index=n;
2717            /*-------------------------------------
2718            if index<n at end of for-loop, then this
2719            means that index and j form a base pair,
2720            and this is noted by numBasePairs[j][i]=index.
2721            if index=n at end of for-loop, then this
2722            means that j is not base paired.
2723            -------------------------------------*/
2724
2725           if ( numBasePairs[i][j-1]>max ){
2726              max = numBasePairs[i][j-1];
2727              index = n;
2728              // j not basepaired with some k such that i<k<j
2729           }
2730
2731           val = basePair(S[i],S[j]) + numBasePairs[i+1][j-1];
2732           if ( j-i<= THRESHOLD && val > max ){
2733              max = val;
2734              index=i;
2735           }
2736           for(k=i; k<=j-THRESHOLD; k++){
2737              val = basePair(S[k],S[j]) + numBasePairs[i][k-1]
2738                    + numBasePairs[k+1][j-1];
2739              if (val > max) {
2740                 max = val;
2741                 index=k;
2742              }
2743           }
2744           numBasePairs[i][j]=max;
2745           if (index<n)
2746              numBasePairs[j][i]=index;
2747           else
2748              numBasePairs[j][i]=-1;
2749        }
2750      }
2751    }
2752
2753 }
2754
2755
2756
2757
2758 static int backtrack(int i, int j, int **numBasePairs,char *S, char *paren, int THRESHOLD)
2759 {
2760   int k;
2761
2762    k = numBasePairs[j][i];
2763    if(k != -1)
2764      {
2765      paren[k] = '(';
2766      paren[j] = ')';
2767      if( THRESHOLD <= (j-1)-(k+1) )
2768        backtrack(k+1,j-1,numBasePairs,S,paren, THRESHOLD);
2769      if (THRESHOLD <= k-1-i  )
2770        backtrack(i,k-1,numBasePairs,S,paren, THRESHOLD);
2771      }
2772    else{
2773      if( THRESHOLD <= j-1-i )
2774        {
2775          backtrack(i,j-1,numBasePairs,S,paren, THRESHOLD);
2776        }
2777      else
2778        return 0;
2779    }
2780    return 0;}
2781
2782 int count;
2783 char * rna_struc2rna_lib ( char *seq_name, char *seq, char *name)
2784 {
2785   FILE *fp;
2786   char *st;
2787
2788
2789   st=nussinov (seq, 2);
2790   if ( name==NULL)name=vtmpnam(NULL);
2791   fp=vfopen ( name, "w");
2792   fprintf (fp, "! TC_LIB_FORMAT_01\n");
2793   fprintf (fp, "1\n%s %d %s\n", seq_name, (int)strlen (seq), seq);
2794   fprintf (fp, "#1 1\n");
2795   display_rna_ss (0, seq, st, fp);
2796   fprintf ( fp, "! SEQ_1_TO_N\n");
2797   vfclose (fp);
2798   vfree (st);
2799   //printf_system ( "cp %s test", name);
2800   return name;
2801 }
2802 int display_rna_ss ( int n, char *seq, char *st, FILE *fp)
2803 {
2804   char p;
2805   char string[100];
2806   static int thread;
2807
2808   while ((p=st[n])!='\0')
2809     {
2810       if ( p=='(')
2811         {
2812           thread=count++;
2813           sprintf (string, "%d",n+1);
2814           n=display_rna_ss (n+1, seq, st, fp);
2815           fprintf (fp, "%s %d 100\n", string, n+1);
2816         }
2817       else if (p=='.');
2818       else if (p==')')
2819         {
2820           return n;
2821         }
2822       n++;
2823     }
2824   return n;
2825 }
2826 /******************************COPYRIGHT NOTICE*******************************/
2827 /*© Centro de Regulacio Genomica */
2828 /*and */
2829 /*Cedric Notredame */
2830 /*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
2831 /*All rights reserved.*/
2832 /*This file is part of T-COFFEE.*/
2833 /**/
2834 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
2835 /*    it under the terms of the GNU General Public License as published by*/
2836 /*    the Free Software Foundation; either version 2 of the License, or*/
2837 /*    (at your option) any later version.*/
2838 /**/
2839 /*    T-COFFEE is distributed in the hope that it will be useful,*/
2840 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
2841 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
2842 /*    GNU General Public License for more details.*/
2843 /**/
2844 /*    You should have received a copy of the GNU General Public License*/
2845 /*    along with Foobar; if not, write to the Free Software*/
2846 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
2847 /*...............................................                                                                                      |*/
2848 /*  If you need some more information*/
2849 /*  cedric.notredame@europe.com*/
2850 /*...............................................                                                                                                                                     |*/
2851 /**/
2852 /**/
2853 /*      */
2854 /******************************COPYRIGHT NOTICE*******************************/