binaries/src/clustalo/src/clustal-omega.c

   1 /* -*- mode: c; tab-width: 4; c-basic-offset: 4;  indent-tabs-mode: nil -*- */
   2
   3 /*********************************************************************
   4  * Clustal Omega - Multiple sequence alignment
   5  *
   6  * Copyright (C) 2010 University College Dublin
   7  *
   8  * Clustal-Omega is free software; you can redistribute it and/or
   9  * modify it under the terms of the GNU General Public License as
  10  * published by the Free Software Foundation; either version 2 of the
  11  * License, or (at your option) any later version.
  12  *
  13  * This file is part of Clustal-Omega.
  14  *
  15  ********************************************************************/
  16
  17 /*
  18  *  RCS $Id: clustal-omega.c 304 2016-06-13 13:39:13Z fabian $
  19  */
  20
  21 #ifdef HAVE_CONFIG_H
  22 #include "config.h"
  23 #endif
  24
  25 #include <assert.h>
  26
  27 #include "clustal-omega.h"
  28 #include "hhalign/general.h"
  29 #include "clustal/hhalign_wrapper.h"
  30
  31 /* The following comment block contains the frontpage/mainpage of the doxygen
  32  *  documentation. Please add some more info. FIXME add more
  33  */
  34
  35 /**
  36  *
  37  * @mainpage Clustal-Omega Documentation
  38  *
  39  * @section intro_sec Introduction
  40  *
  41  * For more information see http://www.clustal.org/
  42  *
  43  * @section api_section API
  44  *
  45  * @subsection example_prog_subsection An Example Program
  46  *
  47  * To use libclustalo you will have to include the clustal-omega.h header and
  48  * link against libclustalo. For linking against libclustalo you will have to
  49  * use a C++ compiler, no matter if your program was written in C or C++. See
  50  * below (\ref pkgconfig_subsubsec)) on how to figure out compiler flags with
  51  * pkg-config.
  52  *
  53  * Assuming Clustal Omega was installed in system-wide default directory (e.g.
  54  * /usr), first compile (don't link yet) your source (for example code see
  55  * section \ref example_src_subsubsec) and then link against libclustalo:
  56  *
  57  * @code
  58  * $ gcc -c -ansi -Wall clustalo-api-test.c
  59  * $ g++  -ansi -Wall -o clustalo-api-test clustalo-api-test.o  -lclustalo
  60  * @endcode
  61  *
  62  * Voila! Now you have your own alignment program based on Clustal Omega which
  63  * can be run with
  64  *
  65  * @code
  66  * $ ./clustalo-api-test <your-sequence-input>
  67  * @endcode
  68  *
  69  * It's best to use the same compiler that you used for compiling libclustal.
  70  * If libclustal was compiled with OpenMP support, you will have to use OpenMP
  71  * flags for you program as well.
  72  *
  73  *
  74  * @subsubsection pkgconfig_subsubsec Using pkg-config / Figuring out compiler flags
  75  *
  76  * Clustal Omega comes with support for <a
  77  * href="http://pkg-config.freedesktop.org">pkg-config</a>, which means you
  78  * can run
  79  *
  80  * @code
  81  * $ pkg-config --cflags --libs clustalo
  82  * @endcode
  83  *
  84  * to figure out cflags and library flags needed to compile and link against
  85  * libclustalo. This is especially handy if Clustal Omega was installed to a
  86  * non-standard directory.
  87  *
  88  * You might have to change PKG_CONFIG_PATH. For example, if you used the prefix $HOME/local/ for
  89  * installation then you will first need to set PKG_CONFIG_PATH:
  90  *
  91  * @code
  92  * $ export PKG_CONFIG_PATH=$HOME/local/lib/pkgconfig
  93  * $ pkg-config --cflags --libs clustalo
  94  * @endcode
  95  *
  96  *
  97  * To compile your source use as above but this time using proper flags:
  98  *
  99  * @code
 100  * $ export PKG_CONFIG_PATH=$HOME/local/lib/pkgconfig
 101  * $ gcc -c -ansi -Wall $(pkg-config --cflags clustalo) clustalo-api-test.c
 102  * $ g++  -ansi -Wall -o clustalo-api-test $(pkg-config --libs clustalo) clustalo-api-test.o
 103  * @endcode
 104  *
 105  *
 106  * @subsubsection example_src_subsubsec Example Source Code
 107  *
 108  * @include "clustalo-api-test.c"
 109  *
 110  *
 111  */
 112
 113
 114 /* the following are temporary flags while the code is still under construction;
 115    had problems internalising hhmake, so as temporary crutch
 116    write alignment to file and get external hmmer/hhmake via system call
 117    to read alignment and convert into HMM
 118    All this will go, once hhmake is properly internalised */
 119 #define INDIRECT_HMM 0 /* temp flag: (1) write aln to file, use system(hmmer/hhmake), (0) internal hhmake */
 120 #define USEHMMER 1 /* temp flag: use system(hmmer) to build HMM */
 121 #define USEHHMAKE (!USEHMMER) /* temp flag: use system(hhmake) to build HMM */
 122
 123
 124 /* shuffle order of input sequences */
 125 #define SHUFFLE_INPUT_SEQ_ORDER 0
 126
 127 /* sort input sequences by length */
 128 #define SORT_INPUT_SEQS 0
 129
 130
 131 int iNumberOfThreads;
 132
 133 /* broken, unused and lonely */
 134 static const int ITERATION_SCORE_IMPROVEMENT_THRESHOLD = 0.01;
 135
 136
 137 /**
 138  * @brief Print Long version information to pre-allocated char.
 139  *
 140  * @note short version
 141  * information is equivalent to PACKAGE_VERSION
 142  *
 143  * @param[out] pcStr
 144  * char pointer to write to preallocated to hold iSize chars.
 145  * @param[in] iSize
 146  * size of pcStr
 147  */
 148 void
 149 PrintLongVersion(char *pcStr, int iSize)
 150 {
 151     snprintf(pcStr, iSize, "version %s; code-name '%s'; build date %s",
 152              PACKAGE_VERSION, PACKAGE_CODENAME, __DATE__);
 153 }
 154 /* end of PrintLongVersion() */
 155
 156
 157
 158 /**
 159  * @brief free aln opts members
 160  *
 161  */
 162 void
 163 FreeAlnOpts(opts_t *prAlnOpts) {
 164     if (NULL != prAlnOpts->pcGuidetreeInfile) {
 165         CKFREE(prAlnOpts->pcGuidetreeInfile);
 166     }
 167     if (NULL != prAlnOpts->pcGuidetreeOutfile) {
 168         CKFREE(prAlnOpts->pcGuidetreeOutfile);
 169     }
 170     if (NULL  != prAlnOpts->pcDistmatOutfile) {
 171         CKFREE(prAlnOpts->pcDistmatOutfile);
 172     }
 173     if (NULL != prAlnOpts->pcDistmatInfile) {
 174         CKFREE(prAlnOpts->pcDistmatInfile);
 175     }
 176 }
 177 /* end of FreeAlnOpts() */
 178
 179
 180
 181 /**
 182  * @brief Sets members of given user opts struct to default values
 183  *
 184  * @param[out] prOpts
 185  * User opt struct to initialise
 186  *
 187  */
 188 void
 189 SetDefaultAlnOpts(opts_t *prOpts) {
 190     prOpts->bAutoOptions = FALSE;
 191
 192     prOpts->pcDistmatInfile = NULL;
 193     prOpts->pcDistmatOutfile = NULL;
 194
 195     prOpts->iClustersizes = 100; /* FS, r274 -> */
 196     prOpts->iTransitivity = 0; /* FS, r290 -> */
 197     prOpts->pcClustfile = NULL; /* FS, r274 -> */
 198     prOpts->pcPosteriorfile = NULL; /* FS, r288 -> */
 199     prOpts->iClusteringType = CLUSTERING_UPGMA;
 200     prOpts->iPairDistType = PAIRDIST_KTUPLE;
 201     prOpts->bUseMbed = TRUE; /* FS, r250 -> */
 202     prOpts->bUseMbedForIteration = TRUE; /* FS, r250 -> */
 203     prOpts->bPileup = FALSE; /* FS, r288 -> */
 204     prOpts->pcGuidetreeOutfile = NULL;
 205     prOpts->pcGuidetreeInfile = NULL;
 206     prOpts->bPercID = FALSE;
 207
 208     prOpts->ppcHMMInput = NULL;
 209     prOpts->iHMMInputFiles = 0;
 210     prOpts->pcHMMBatch = NULL; /* FS, r291 -> */
 211
 212     prOpts->iNumIterations = 0;
 213     prOpts->bIterationsAuto = FALSE;
 214     prOpts->iMaxGuidetreeIterations = INT_MAX;
 215     prOpts->iMaxHMMIterations = INT_MAX;
 216
 217     SetDefaultHhalignPara(& prOpts->rHhalignPara);
 218  }
 219 /* end of SetDefaultAlnOpts() */
 220
 221
 222
 223 /**
 224  * @brief Check logic of parsed user options. Will exit (call Log(&rLog, LOG_FATAL, ))
 225  * on Fatal logic error
 226  *
 227  * @param[in] prOpts
 228  * Already parsed user options
 229  *
 230  */
 231 void
 232 AlnOptsLogicCheck(opts_t *prOpts)
 233 {
 234     /* guide-tree & distmat
 235      *
 236      */
 237     if (prOpts->pcDistmatInfile && prOpts->pcGuidetreeInfile) {
 238         Log(&rLog, LOG_FATAL, "Read distances *and* guide-tree from file doesn't make sense.");
 239     }
 240
 241     if (prOpts->pcDistmatOutfile && prOpts->pcGuidetreeInfile) {
 242         Log(&rLog, LOG_FATAL, "Won't be able to save distances to file, because I got a guide-tree as input.");
 243     }
 244
 245     /* combination of options that don't make sense when not iterating
 246      */
 247     if (prOpts->iNumIterations==0 && prOpts->bIterationsAuto != TRUE)  {
 248
 249         if (prOpts->pcGuidetreeInfile && prOpts->pcGuidetreeOutfile) {
 250             Log(&rLog, LOG_FATAL, "Got a guide-tree as input and output which doesn't make sense when not iterating.");
 251         }
 252         /*
 253           if (prOpts->pcGuidetreeInfile && prOpts->bUseMbed > 0) {
 254           Log(&rLog, LOG_FATAL, "Got a guide-tree as input and was requested to cluster with mBed, which doesn't make sense when not iterating.");
 255           }
 256         */
 257         /*
 258           AW: bUseMbedForIteration default since at least R252
 259           if (prOpts->bUseMbedForIteration > 0) {
 260             Log(&rLog, LOG_FATAL, "No iteration requested, but mbed for iteration was set. Paranoia exit.");
 261           }
 262         */
 263     }
 264
 265     if (prOpts->rHhalignPara.iMacRamMB < 512) {
 266         Log(&rLog, LOG_WARN, "Memory for MAC Algorithm quite low, Viterbi Algorithm may be triggered.");
 267         if (prOpts->rHhalignPara.iMacRamMB < 1) {
 268             Log(&rLog, LOG_WARN, "Viterbi Algorithm always turned on, increase MAC-RAM to turn on MAC.");
 269         }
 270     }
 271
 272     return;
 273 }
 274 /* end of AlnOptsLogicCheck() */
 275
 276
 277 /**
 278  * @brief FIXME doc
 279  */
 280 void
 281 PrintAlnOpts(FILE *prFile, opts_t *prOpts)
 282 {
 283     int iAux;
 284
 285
 286     /* keep in same order as struct */
 287     fprintf(prFile, "option: auto-options = %d\n", prOpts->bAutoOptions);
 288     fprintf(prFile, "option: distmat-infile = %s\n",
 289             NULL != prOpts->pcDistmatInfile? prOpts->pcDistmatInfile: "(null)");
 290         fprintf(prFile, "option: distmat-outfile = %s\n",
 291             NULL != prOpts->pcDistmatOutfile? prOpts->pcDistmatOutfile: "(null)");
 292         fprintf(prFile, "option: clustering-type = %d\n", prOpts->iClusteringType);
 293         fprintf(prFile, "option: pair-dist-type = %d\n", prOpts->iPairDistType);
 294         fprintf(prFile, "option: use-mbed = %d\n", prOpts->bUseMbed);
 295         fprintf(prFile, "option: use-mbed-for-iteration = %d\n", prOpts->bUseMbedForIteration);
 296         fprintf(prFile, "option: pile-up = %d\n", prOpts->bPileup);
 297         fprintf(prFile, "option: guidetree-outfile = %s\n",
 298             NULL != prOpts->pcGuidetreeOutfile? prOpts->pcGuidetreeOutfile: "(null)");
 299         fprintf(prFile, "option: guidetree-infile = %s\n",
 300             NULL != prOpts->pcGuidetreeInfile? prOpts->pcGuidetreeInfile: "(null)");
 301     for (iAux=0; iAux<prOpts->iHMMInputFiles; iAux++) {
 302         fprintf(prFile, "option: hmm-input no %d = %s\n", iAux, prOpts->ppcHMMInput[iAux]);
 303     }
 304         fprintf(prFile, "option: hmm-input-files = %d\n", prOpts->iHMMInputFiles);
 305         fprintf(prFile, "option: num-iterations = %d\n", prOpts->iNumIterations);
 306         fprintf(prFile, "option: iterations-auto = %d\n", prOpts->bIterationsAuto);
 307         fprintf(prFile, "option: max-hmm-iterations = %d\n", prOpts->iMaxHMMIterations);
 308         fprintf(prFile, "option: max-guidetree-iterations = %d\n", prOpts->iMaxGuidetreeIterations);
 309     fprintf(prFile, "option: iMacRamMB = %d\n", prOpts->rHhalignPara.iMacRamMB);
 310     fprintf(prFile, "option: percent-id = %d\n", prOpts->bPercID);
 311     fprintf(prFile, "option: use-kimura = %d\n", prOpts->bUseKimura);
 312     fprintf(prFile, "option: clustering-out = %s\n", prOpts->pcClustfile);
 313     fprintf(prFile, "option: posterior-out = %s\n", prOpts->pcPosteriorfile);
 314
 315 }
 316 /* end of PrintAlnOpts() */
 317
 318
 319
 320 /**
 321  * @brief Returns major version of HMMER. Whichever hmmbuild version
 322  * is found first in your PATH will be used
 323  *
 324  * @return -1 on error, major hmmer version otherwise
 325  *
 326  */
 327 int
 328 HmmerVersion()
 329 {
 330     char zcHmmerTestCall[] = "hmmbuild -h";
 331     FILE *fp = NULL;
 332     int iMajorVersion = 0;
 333     char zcLine[16384];
 334
 335     if (NULL == (fp = popen(zcHmmerTestCall, "r"))) {
 336         Log(&rLog, LOG_ERROR, "Couldn't exec %s", zcHmmerTestCall);
 337         return -1;
 338     }
 339     while (fgets(zcLine, sizeof(zcLine), fp)) {
 340         char *pcLocate;
 341         if ((pcLocate = strstr(zcLine, "HMMER "))) {
 342             iMajorVersion = atoi(&pcLocate[6]);
 343             break;
 344         }
 345     }
 346     pclose(fp);
 347
 348     return iMajorVersion;
 349 }
 350 /* end of HmmerVersion() */
 351
 352
 353
 354 /**
 355  * @brief Create a HHM file from aligned sequences
 356  *
 357  * @warning Should be eliminated in the future
 358  * as building routine should not create intermediate files
 359  *
 360  * @param[in] prMSeq
 361  * Aligned mseq_t
 362  * @param[in] pcHMMOut
 363  * HMM output file name
 364  *
 365  * @return Non-zero on error
 366  *
 367  */
 368 int
 369 AlnToHHMFile(mseq_t *prMSeq, char *pcHMMOut)
 370 {
 371     char *tmp_aln = NULL;
 372     int retcode = OK;
 373
 374     assert(NULL!=prMSeq);
 375     assert(NULL!=pcHMMOut);
 376
 377     if (FALSE == prMSeq->aligned) {
 378         Log(&rLog, LOG_ERROR, "Sequences need to be aligned to create an HMM");
 379         return FAILURE;
 380     }
 381
 382     /* Convert alignment to a2m, and call hhmake
 383      *
 384      * can't be static templates, or mktemp fails (at least on os x
 385      * (with a bus error))
 386      *
 387      * gcc says we should use mkstemp to avoid race conditions,
 388      * but that returns a file descriptor, which is of no use to
 389      * us
 390      */
 391     /* NOTE: the following won't work on windows: missing /tmp/ */
 392     tmp_aln = CkStrdup("/tmp/clustalo_tmpaln_XXXXXX");
 393     if (NULL == mktemp(tmp_aln)) {
 394         Log(&rLog, LOG_ERROR, "Could not create temporary alignment filename");
 395         retcode = FAILURE;
 396         goto cleanup_and_return;
 397     }
 398 #define LINE_WRAP 60
 399     if (WriteAlignment(prMSeq, tmp_aln, MSAFILE_A2M, LINE_WRAP, FALSE)) {
 400         Log(&rLog, LOG_ERROR, "Could not save alignment to %s", tmp_aln);
 401         retcode = FAILURE;
 402         goto cleanup_and_return;
 403     }
 404
 405     if (HHMake_Wrapper(tmp_aln, pcHMMOut)){
 406         Log(&rLog, LOG_ERROR, "Could not convert alignment %s into HHM", tmp_aln);
 407         retcode = FAILURE;
 408         goto cleanup_and_return;
 409     }
 410
 411
 412  cleanup_and_return:
 413
 414     if (NULL != tmp_aln) {
 415         if (FileExists(tmp_aln)) {
 416             if (remove(tmp_aln)) {
 417                 Log(&rLog, LOG_WARN, "Removing %s failed. Continuing anyway", tmp_aln);
 418             }
 419         }
 420         CKFREE(tmp_aln);
 421     }
 422
 423     return retcode;
 424
 425 } /* end of AlnToHHMFile() */
 426
 427
 428
 429 /**
 430  * @brief Create a HMM file from aligned sequences
 431  *
 432  * @warning Should be replaced in the future by some internal HMM
 433  * building routine that does not call external programs
 434  *
 435  * @param[in] prMSeq
 436  * Aligned mseq_t
 437  * @param[in] pcHMMOut
 438  * HMM output file name
 439  *
 440  * @return Non-zero on error
 441  *
 442
 443  */
 444 int
 445 AlnToHMMFile(mseq_t *prMSeq, const char *pcHMMOut)
 446 {
 447     char *tmp_aln = NULL;
 448     char *tmp_hmm = NULL; /* only needed for hmmer3 to hmmer2 conversion */
 449     char cmdbuf[16384];
 450     int iHmmerVersion = 0;
 451     int retcode = OK;
 452
 453     assert(NULL!=prMSeq);
 454     assert(NULL!=pcHMMOut);
 455
 456     if (FALSE == prMSeq->aligned) {
 457         Log(&rLog, LOG_ERROR, "Sequences need to be aligned to create an HMM");
 458         return FAILURE;
 459     }
 460
 461     iHmmerVersion = HmmerVersion();
 462     if (2 != iHmmerVersion && 3 != iHmmerVersion) {
 463         Log(&rLog, LOG_ERROR, "Could not find suitable HMMER binaries");
 464         return FAILURE;
 465     }
 466
 467     /* Convert alignment to stockholm, call hmmbuild and then
 468      * either hmmconvert (hmmer3) or hmmcalibrate (hmmer2)
 469      *
 470      * can't be static templates, or mktemp fails (at least on os x
 471      * (with a bus error))
 472      *
 473      * gcc says we should use mkstemp to avoid race conditions,
 474      * but that returns a file descriptor, which is of no use to
 475      * us
 476      */
 477     /* NOTE: the following won't work on windows: missing /tmp/ */
 478     tmp_aln = CkStrdup("/tmp/clustalo_tmpaln_XXXXXX");
 479     if (NULL == mktemp(tmp_aln)) {
 480         Log(&rLog, LOG_ERROR, "Could not create temporary alignment filename");
 481         retcode = FAILURE;
 482         goto cleanup_and_return;
 483     }
 484 #define LINE_WRAP 60
 485     if (WriteAlignment(prMSeq, tmp_aln, MSAFILE_STOCKHOLM, LINE_WRAP, FALSE)) {
 486         Log(&rLog, LOG_ERROR, "Could not save alignment to %s", tmp_aln);
 487         retcode = FAILURE;
 488         goto cleanup_and_return;
 489     }
 490
 491     if (2 == iHmmerVersion) {
 492         sprintf(cmdbuf, "hmmbuild %s %s >/dev/null && hmmcalibrate %s >/dev/null",
 493                 pcHMMOut, tmp_aln, pcHMMOut);
 494         if (system(cmdbuf)) {
 495             Log(&rLog, LOG_ERROR, "Command '%s' failed", cmdbuf);
 496             retcode = FAILURE;
 497             goto cleanup_and_return;
 498         }
 499     } else if (3 == iHmmerVersion) {
 500         /* NOTE: the following won't work on windows: missing /tmp/ */
 501         tmp_hmm = CkStrdup("/tmp/clustalo_tmphmm2_XXXXXX");
 502         if (NULL == mktemp(tmp_hmm)) {
 503             Log(&rLog, LOG_ERROR, "Could not create temporary hmm filename");
 504             retcode = FAILURE;
 505             goto cleanup_and_return;
 506         }
 507         sprintf(cmdbuf, "hmmbuild %s %s >/dev/null && hmmconvert -2 %s > %s",
 508                 tmp_hmm, tmp_aln, tmp_hmm, pcHMMOut);
 509         if (system(cmdbuf)) {
 510             Log(&rLog, LOG_ERROR, "Command '%s' failed", cmdbuf);
 511             retcode = FAILURE;
 512             goto cleanup_and_return;
 513         }
 514     } else {
 515         CKFREE(tmp_aln);
 516         Log(&rLog, LOG_FATAL, "Internal error: Unknown Hmmer version %d", iHmmerVersion);
 517     }
 518
 519
 520  cleanup_and_return:
 521
 522     if (NULL != tmp_aln) {
 523         if (FileExists(tmp_aln)) {
 524             if (remove(tmp_aln)) {
 525                 Log(&rLog, LOG_WARN, "Removing %s failed. Continuing anyway", tmp_aln);
 526             }
 527         }
 528         CKFREE(tmp_aln);
 529     }
 530     if (NULL != tmp_hmm) {
 531         if (FileExists(tmp_hmm)) {
 532             if (remove(tmp_hmm)) {
 533                 Log(&rLog, LOG_WARN, "Removing %s failed. Continuing anyway", tmp_hmm);
 534             }
 535         }
 536         CKFREE(tmp_hmm);
 537     }
 538
 539     return retcode;
 540 }
 541 /* end of AlnToHMMFile() */
 542
 543
 544
 545 /**
 546  * @brief Convert a multiple sequence structure into a HMM
 547  *
 548  * @param[out] prHMM
 549  * Pointer to preallocted HMM which will be set here
 550  * @param[in] prMSeq
 551  * Pointer to an alignment
 552  *
 553  * @return 0 on error, non-0 otherwise
 554  *
 555  * @see AlnToHMMFile()
 556  *
 557  */
 558 int
 559 AlnToHMM(hmm_light *prHMM, mseq_t *prMSeq)
 560 {
 561     char *pcHMM; /* temp hmm file */
 562
 563     Log(&rLog, LOG_INFO,
 564         "Using HMMER version %d to calculate a new HMM.",
 565          HmmerVersion());
 566     /* FIXME replace all this with internal HMM computation (HHmake) */
 567
 568     /**
 569      * @warning the following probably won't work on windows: missing
 570      * /tmp/. Should be ok on Cygwin though
 571      */
 572     pcHMM = CkStrdup("/tmp/clustalo-hmm-iter_XXXXXX");
 573     if (NULL == mktemp(pcHMM)) {
 574         Log(&rLog, LOG_ERROR, "Could not create temporary hmm filename");
 575         CKFREE(pcHMM);
 576         return FAILURE;
 577     }
 578
 579     /* Create a HMM representing the current alignment
 580      */
 581 #if USEHMMER
 582     if (AlnToHMMFile(prMSeq, pcHMM)) {
 583         Log(&rLog, LOG_ERROR, "AlnToHMMFile() on %s failed.", pcHMM);
 584         CKFREE(pcHMM);
 585         return FAILURE;
 586     }
 587 #elif USEHHMAKE
 588     if (AlnToHHMFile(prMSeq, pcHMM)) {
 589         Log(&rLog, LOG_ERROR, "AlnToHHMFile() on %s failed.", pcHMM);
 590         CKFREE(pcHMM);
 591         return FAILURE;
 592     }
 593     /* Log(&rLog, LOG_FATAL, "Method to create HHM (HMM using hhmake) not installed yet"); */
 594 #else
 595     Log(&rLog, LOG_FATAL, "Unknown method to create temporary HMM");
 596 #endif
 597
 598     /* Read HMM information
 599      */
 600     if (OK != readHMMWrapper(prHMM, pcHMM)){
 601         Log(&rLog, LOG_ERROR, "Processing of HMM file %s failed", pcHMM);
 602         CKFREE(pcHMM);
 603         return FAILURE;
 604     }
 605
 606     if (remove(pcHMM)) {
 607         Log(&rLog, LOG_WARN, "Removing %s failed. Continuing anyway", pcHMM);
 608     }
 609     CKFREE(pcHMM);
 610
 611     return OK;
 612 }
 613 /* end of AlnToHMM() */
 614
 615
 616
 617 /**
 618  * @brief FIXME
 619  *
 620  */
 621 void
 622 InitClustalOmega(int iNumThreadsRequested)
 623 {
 624
 625 #ifdef HAVE_OPENMP
 626     iNumberOfThreads = iNumThreadsRequested;
 627     omp_set_num_threads(iNumberOfThreads);
 628 #else
 629     if (iNumThreadsRequested>1) {
 630         Log(&rLog, LOG_FATAL, "Cannot change number of threads to %d. %s was build without OpenMP support.",
 631               iNumThreadsRequested, PACKAGE_NAME);
 632     }
 633     iNumberOfThreads = 1; /* need to set this, even if build without support */
 634 #endif
 635
 636     Log(&rLog, LOG_INFO, "Using %d threads",
 637          iNumberOfThreads);
 638
 639 }
 640 /* end of InitClustalOmega() */
 641
 642
 643
 644 /**
 645  * @brief Defines an alignment order, which adds sequences
 646  * sequentially, i.e. one at a time starting with seq 1 & 2
 647  *
 648  * @param[out] piOrderLR_p
 649  * order in which nodes/profiles are to be merged/aligned
 650  * @param[in] iNumSeq
 651  * Number of sequences
 652  *
 653  * @see TraverseTree()
 654  *
 655  */
 656 void
 657 SequentialAlignmentOrder(int **piOrderLR_p, int iNumSeq)
 658 {
 659     unsigned int uNodes = iNumSeq*2-1;
 660     unsigned int uNodeCounter = 0;
 661     unsigned int uSeqCounter = 0;
 662
 663     Log(&rLog, LOG_FATAL, "FIXME: Untested...");
 664
 665     (*piOrderLR_p) = (int *)CKCALLOC(DIFF_NODE * uNodes, sizeof(int));
 666     /* loop over merge nodes, which have per definition even indices
 667      * and set up children which have odd indices
 668      */
 669     uSeqCounter = 0;
 670     for (uNodeCounter=iNumSeq; uNodeCounter<uNodes; uNodeCounter+=1) {
 671          unsigned int uLeftChildNodeIndex = uNodeCounter-1;
 672          unsigned int uRightChildNodeIndex = uNodeCounter-iNumSeq+1;
 673          unsigned int uParentNodeIndex = uNodeCounter+1;
 674
 675          /* merge node setup */
 676         (*piOrderLR_p)[DIFF_NODE*uNodeCounter+LEFT_NODE] = uLeftChildNodeIndex;
 677         (*piOrderLR_p)[DIFF_NODE*uNodeCounter+RGHT_NODE] = uRightChildNodeIndex;
 678         (*piOrderLR_p)[DIFF_NODE*uNodeCounter+PRNT_NODE] = uParentNodeIndex;
 679         /* only setup left child if at first merge node, all other left childs
 680          * should be merge nodes that are already set up. also correct
 681          * left node number here.
 682          */
 683         if (uNodeCounter==iNumSeq) {
 684             (*piOrderLR_p)[DIFF_NODE*uNodeCounter+LEFT_NODE] = 0;
 685
 686             (*piOrderLR_p)[0+LEFT_NODE] = 0;
 687             (*piOrderLR_p)[0+RGHT_NODE] = 0;
 688             (*piOrderLR_p)[0+PRNT_NODE] = uNodeCounter;
 689             uSeqCounter++;
 690
 691             Log(&rLog, LOG_FORCED_DEBUG, "Set up first leaf with node counter %d: left=%d right=%d parent=%d",
 692                       0,
 693                       (*piOrderLR_p)[DIFF_NODE*uLeftChildNodeIndex+LEFT_NODE],
 694                       (*piOrderLR_p)[DIFF_NODE*uLeftChildNodeIndex+RGHT_NODE],
 695                       (*piOrderLR_p)[DIFF_NODE*uLeftChildNodeIndex+PRNT_NODE]);
 696         }
 697         Log(&rLog, LOG_FORCED_DEBUG, "Set up merge node with node counter %d: left=%d right=%d parent=%d",
 698                   uNodeCounter, (*piOrderLR_p)[DIFF_NODE*uNodeCounter+LEFT_NODE],
 699                   (*piOrderLR_p)[DIFF_NODE*uNodeCounter+RGHT_NODE],
 700                   (*piOrderLR_p)[DIFF_NODE*uNodeCounter+PRNT_NODE]);
 701
 702         /* right child */
 703         (*piOrderLR_p)[DIFF_NODE*uRightChildNodeIndex+LEFT_NODE] = uSeqCounter;
 704         (*piOrderLR_p)[DIFF_NODE*uRightChildNodeIndex+RGHT_NODE] = uSeqCounter;
 705         (*piOrderLR_p)[DIFF_NODE*uRightChildNodeIndex+PRNT_NODE] = uNodeCounter;
 706         uSeqCounter++;
 707
 708         Log(&rLog, LOG_FORCED_DEBUG, "Set up leaf with node counter %d: left=%d right=%d parent=%d",
 709                   uRightChildNodeIndex, (*piOrderLR_p)[DIFF_NODE*uRightChildNodeIndex+LEFT_NODE],
 710                   (*piOrderLR_p)[DIFF_NODE*uRightChildNodeIndex+RGHT_NODE],
 711                   (*piOrderLR_p)[DIFF_NODE*uRightChildNodeIndex+PRNT_NODE]);
 712     }
 713 }
 714 /* end of SequentialAlignmentOrder() */
 715
 716
 717
 718 /**
 719  * @brief Defines the alignment order by calculating a guide tree. In
 720  * a first-step pairwise distances will be calculated (or read from a
 721  * file). In a second step those distances will be clustered and a
 722  * guide-tree created. Steps 1 and 2 will be skipped if a guide-tree
 723  * file was given, in which case the guide-tree will be just read from
 724  * the file.
 725  *
 726  * @param[out] piOrderLR_p
 727  * order in which nodes/profiles are to be merged/aligned
 728  * @param[out] pdSeqWeights_p
 729  * Sequence weights
 730  * @param[out] pdSeqWeights_p
 731  * Sequence weights
 732  * @param[in] prMSeq
 733  * The sequences from which the alignment order is to be calculated
 734  * @param[in] iPairDistType
 735  * Method of pairwise distance comparison
 736  * @param[in] pcDistmatInfile
 737  * If not NULL distances will be read from this file instead of being
 738  * calculated
 739  * @param[in] pcDistmatOutfile
 740  * If not NULL computed pairwise distances will be written to this file
 741  * @param[in] iClusteringType
 742  * Clustering method to be used to cluster the pairwise distances
 743  * @param[in] pcGuidetreeInfile
 744  * If not NULL guidetree will be read from this file. Skips pairwise
 745  * distance and guidetree computation
 746  * @param[in] pcGuidetreeOutfile
 747  * If not NULL computed guidetree will be written to this file
 748  * @param[in] bUseMbed
 749  * If TRUE, fast mBed guidetree computation will be employed
 750  *
 751  * @return Non-zero on error
 752  *
 753  */
 754 int
 755 AlignmentOrder(int **piOrderLR_p, double **pdSeqWeights_p, mseq_t *prMSeq,
 756                int iPairDistType, char *pcDistmatInfile, char *pcDistmatOutfile,
 757                int iClusteringType, int iClustersizes,
 758                char *pcGuidetreeInfile, char *pcGuidetreeOutfile, char *pcClusterFile,
 759                bool bUseMbed, bool bPercID)
 760 {
 761     /* pairwise distance matrix (tmat in 1.83) */
 762     symmatrix_t *distmat = NULL;
 763     /* guide tree */
 764     tree_t *prTree = NULL;
 765     int i = 0;
 766
 767
 768     /* Shortcut for only two sequences: Do not compute k-tuple
 769      * distances. Use the same logic as in TraverseTree() to setup
 770      * piOrderLR_p. Changes there will have to be reflected here as
 771      * well. */
 772     if (2==prMSeq->nseqs) {
 773         Log(&rLog, LOG_VERBOSE,
 774             "Have only two sequences: No need to compute pairwise score and compute a tree.");
 775
 776         if (NULL != pcDistmatOutfile){
 777             Log(&rLog, LOG_WARN, "Have only two sequences: Will not calculate/print distance matrix.");
 778         }
 779
 780         (*piOrderLR_p) = (int*) CKMALLOC(DIFF_NODE * 3 * sizeof(int));
 781         (*piOrderLR_p)[DIFF_NODE*0+LEFT_NODE] = 0;
 782         (*piOrderLR_p)[DIFF_NODE*0+RGHT_NODE] = 0;
 783         (*piOrderLR_p)[DIFF_NODE*0+PRNT_NODE] = 0;
 784
 785         (*piOrderLR_p)[DIFF_NODE*1+LEFT_NODE] = 1;
 786         (*piOrderLR_p)[DIFF_NODE*1+RGHT_NODE] = 1;
 787         (*piOrderLR_p)[DIFF_NODE*1+PRNT_NODE] = 1;
 788
 789         /* root */
 790         (*piOrderLR_p)[DIFF_NODE*2+LEFT_NODE] = 0;
 791         (*piOrderLR_p)[DIFF_NODE*2+RGHT_NODE] = 1;
 792         (*piOrderLR_p)[DIFF_NODE*2+PRNT_NODE] = 2;
 793
 794         /* Same logic as CalcClustalWeights(). Changes there will
 795            have to be reflected here as well. */
 796 #if USE_WEIGHTS
 797          (*pdWeights_p) = (double *) CKMALLOC(uNodeCount * sizeof(double));
 798          (*pdWeights_p)[0] = 0.5;
 799          (*pdWeights_p)[1] = 0.5;
 800 #endif
 801
 802         return OK;
 803     }
 804
 805
 806     /* compute distance & guide tree, alternatively read distances or
 807      * guide tree from file
 808      *
 809      */
 810     if (NULL != pcGuidetreeInfile) {
 811         Log(&rLog, LOG_INFO, "Reading guide-tree from %s", pcGuidetreeInfile);
 812         if (GuideTreeFromFile(&prTree, prMSeq, pcGuidetreeInfile)) {
 813             Log(&rLog, LOG_ERROR, "Reading of guide tree %s failed.", pcGuidetreeInfile);
 814             return FAILURE;
 815         }
 816
 817     } else {
 818
 819         if (bUseMbed) {
 820             if (NULL != pcDistmatInfile) {
 821                 Log(&rLog, LOG_ERROR, "Can't input distance matrix when in mbed mode.");
 822                 return FAILURE;
 823             }
 824             if (Mbed(&prTree, prMSeq, iPairDistType, pcGuidetreeOutfile, iClustersizes, pcClusterFile)) {
 825                 Log(&rLog, LOG_ERROR, "mbed execution failed.");
 826                 return FAILURE;
 827             }
 828             Log(&rLog, LOG_INFO, "Guide-tree computation (mBed) done.");
 829             if (NULL != pcDistmatOutfile) {
 830                 Log(&rLog, LOG_INFO,
 831                     "Ignoring request to write distance matrix (am in mBed mode)");
 832             }
 833         } else {
 834
 835             if (PairDistances(&distmat, prMSeq, iPairDistType, bPercID,
 836                               0, prMSeq->nseqs, 0, prMSeq->nseqs,
 837                               pcDistmatInfile, pcDistmatOutfile)) {
 838                 Log(&rLog, LOG_ERROR, "Couldn't compute pair distances");
 839                 return FAILURE;
 840             }
 841
 842             /* clustering of distances to get guide tree
 843              */
 844             if (CLUSTERING_UPGMA == iClusteringType) {
 845                 char **labels;
 846                 labels = (char**) CKMALLOC(prMSeq->nseqs * sizeof(char*));
 847                 for (i=0; i<prMSeq->nseqs; i++) {
 848                     labels[i] = prMSeq->sqinfo[i].name;
 849                 }
 850
 851                 GuideTreeUpgma(&prTree, labels, distmat, pcGuidetreeOutfile);
 852                 Log(&rLog, LOG_INFO, "Guide-tree computation done.");
 853
 854                 CKFREE(labels);
 855             } else {
 856                 Log(&rLog, LOG_FATAL, "INTERNAL ERROR %s",
 857                       "clustering method should have been checked before");
 858             }
 859         } /* did not use mBed */
 860     } /* had to calculate tree (did not read from file) */
 861
 862
 863 #if USE_WEIGHTS
 864     /* derive sequence weights from tree
 865      *
 866      */
 867     Log(&rLog, LOG_INFO, "Calculating sequence weights");
 868     CalcClustalWeights(pdSeqWeights_p, prTree);
 869     for (i = 0; i < GetLeafCount(prTree); i++) {
 870         Log(&rLog, LOG_VERBOSE,
 871             "Weight for seq no %d: %s = %f",
 872              i, prMSeq->sqinfo[i].name, (*pdSeqWeights_p)[i]);
 873     }
 874 #else
 875     Log(&rLog, LOG_DEBUG, "Not using weights");
 876 #endif
 877
 878
 879     /* define traversing order of tree
 880      *
 881      */
 882     TraverseTree(piOrderLR_p, prTree, prMSeq);
 883     if (rLog.iLogLevelEnabled <= LOG_DEBUG) {
 884         /* FIXME: debug only, FS */
 885         uint uNodeIndex;
 886         FILE *fp = LogGetFP(&rLog, LOG_INFO);
 887         Log(&rLog, LOG_DEBUG, "left/right order after tree traversal");
 888         for (uNodeIndex = 0; uNodeIndex < GetNodeCount(prTree); uNodeIndex++) {
 889             fprintf(fp, "%3d:\t%2d/%2d -> %d\n", i,
 890                    (*piOrderLR_p)[DIFF_NODE*uNodeIndex+LEFT_NODE],
 891                    (*piOrderLR_p)[DIFF_NODE*uNodeIndex+RGHT_NODE],
 892                    (*piOrderLR_p)[DIFF_NODE*uNodeIndex+PRNT_NODE]);
 893         }
 894     }
 895
 896     FreeMuscleTree(prTree);
 897     FreeSymMatrix(&distmat);
 898
 899 #if 0
 900     Log(&rLog, LOG_FATAL, "DEBUG EXIT before leaving %s", __FUNCTION__);
 901 #endif
 902     return OK;
 903 }
 904 /* end of AlignmentOrder() */
 905
 906
 907
 908 /**
 909  * @brief Set some options automatically based on number of sequences. Might
 910  * overwrite some user-set options.
 911  *
 912  * @param[out] prOpts
 913  * Pointer to alignment options structure
 914  * @param[in] iNumSeq
 915  * Number of sequences to align
 916  */
 917 void
 918 SetAutoOptions(opts_t *prOpts, int iNumSeq) {
 919
 920     Log(&rLog, LOG_INFO,
 921         "Setting options automatically based on input sequence characteristics (might overwrite some of your options).");
 922
 923     /* AW: new version of mbed is always good (uses subclusters) */
 924     if (FALSE == prOpts->bUseMbed) {
 925         Log(&rLog, LOG_INFO, "Auto settings: Enabling mBed.");
 926         prOpts->bUseMbed = TRUE;
 927     }
 928
 929     if (iNumSeq >= 1000) {
 930         if (0 != prOpts->iNumIterations) {
 931             Log(&rLog, LOG_INFO, "Auto settings: Disabling iterations.");
 932             prOpts->iNumIterations = 0;
 933         }
 934
 935     } else if (iNumSeq < 1000) {
 936         if (1 != prOpts->iNumIterations) {
 937             Log(&rLog, LOG_INFO, "Auto settings: Setting iteration to 1.");
 938             prOpts->iNumIterations = 1;
 939         }
 940     }
 941 }
 942 /* end of */
 943
 944
 945
 946 /**
 947  * @brief The main alignment function which wraps everything else.
 948  *
 949  * @param[out] prMSeq
 950  * *the* multiple sequences structure
 951  * @param[in] prMSeqProfile
 952  * optional profile to align against
 953  * @param[in] prOpts
 954  * alignment options to use
 955  *
 956  * @return 0 on success, -1 on failure
 957  *
 958  */
 959 int
 960 Align(mseq_t *prMSeq,
 961       mseq_t *prMSeqProfile,
 962       opts_t *prOpts) { /* Note DEVEL 291: at this stage pppcHMMBNames is set but ppiHMMBindex is not */
 963
 964     /* HMM
 965      */
 966     /* structs with pseudocounts etc; one for each HMM infile, i.e.
 967      * index range: 0..iHMMInputFiles */
 968     hmm_light *prHMMs = NULL;
 969
 970     /* MSA order in which nodes/profiles are to be merged/aligned
 971        (order of nodes in guide tree (left/right)*/
 972     int *piOrderLR = NULL;
 973
 974     /* weights per sequence */
 975     double *pdSeqWeights = NULL;
 976
 977     /* Iteration
 978      */
 979     int iIterationCounter = 0;
 980     double dAlnScore;
 981     /* last dAlnScore for iteration */
 982     double dLastAlnScore = -666.666;
 983
 984     /* HMM batch file */
 985     char **ppcHMMbatch = NULL; /* names of unique HMM files */
 986     int iHMMbatch = 0; /* number of unique HMM files */
 987
 988     int i, j; /* aux */
 989
 990     assert(NULL != prMSeq);
 991     if (NULL != prMSeqProfile) {
 992         assert(TRUE == prMSeqProfile->aligned);
 993     }
 994
 995
 996     /* automatic setting of options
 997      *
 998      */
 999     if (prOpts->bAutoOptions) {
1000         SetAutoOptions(prOpts, prMSeq->nseqs);
1001     }
1002
1003
1004 #if SHUFFLE_INPUT_SEQ_ORDER
1005     /*
1006      * shuffle input:  only useful for testing/debugging
1007      */
1008     Log(&rLog, LOG_WARN, "Shuffling input sequences! (Will also change output order)");
1009     ShuffleMSeq(prMSeq);
1010 #endif
1011
1012
1013 #if SORT_INPUT_SEQS
1014     /*
1015      * sort input:
1016      *
1017      * would ensure we *always* (unless we get into the mbed k-means stage)
1018      * get the same answer. usually you don't, because most pairwise alignment
1019      * scores are in theory not symmetric, therefore sequence ordering might
1020      * have an effect on the guide-tree. Sorting by length should get rid of
1021      * this (and takes no time even for 100k seqs). Benchmark results on
1022      * Balibase show almost no difference after sorting.
1023      */
1024     Log(&rLog, LOG_WARN, "Sorting input seq by length! This will also change the output order");
1025     SortMSeqByLength(prMSeq, 'd');
1026
1027 #endif
1028
1029     if (TRUE == prOpts->bPileup){
1030         PileUp(prMSeq, prOpts->rHhalignPara, prOpts->iClustersizes);
1031         return 0;
1032     }
1033
1034
1035     /* Read backgrounds HMMs and store in prHMMs (Devel 291)
1036      *
1037      */
1038     if (NULL != prOpts->pcHMMBatch){
1039         int i, j, k;
1040
1041         for (i = 0; i < prMSeq->nseqs; i++){
1042
1043             if (NULL != prMSeq->pppcHMMBNames[i]){
1044                 for (j = 0; NULL != prMSeq->pppcHMMBNames[i][j]; j++){
1045
1046                     for (k = 0; k < iHMMbatch; k++){
1047                         if (0 == strcmp(ppcHMMbatch[k], prMSeq->pppcHMMBNames[i][j])){
1048                             prMSeq->ppiHMMBindex[i][j] = k;
1049                             break; /* HMM already registered */
1050                         }
1051                     } /* went through HMM batch files already identified */
1052                     if (k == iHMMbatch){
1053                         FILE *pfHMM = NULL;
1054                         if (NULL == (pfHMM = fopen(prMSeq->pppcHMMBNames[i][j], "r"))){
1055                             prMSeq->ppiHMMBindex[i][j] = -1;
1056                             Log(&rLog, LOG_WARN, "Background HMM %s for %s (%d/%d) does not exist",
1057                                 prMSeq->pppcHMMBNames[i][j], prMSeq->sqinfo[i].name, i, j);
1058                         }
1059                         else {
1060                             fclose(pfHMM); pfHMM = NULL;
1061                             ppcHMMbatch = (char **)realloc(ppcHMMbatch, (iHMMbatch+1)*sizeof(char *));
1062                             ppcHMMbatch[iHMMbatch] = strdup(prMSeq->pppcHMMBNames[i][j]);
1063                             prMSeq->ppiHMMBindex[i][j] = k;
1064                             iHMMbatch++;
1065                         }
1066                     }
1067
1068                 } /* j = 0; NULL != prMSeq->pppcHMMBNames[i][j] */
1069             } /* NULL != prMSeq->pppcHMMBNames[i] */
1070             else {
1071                 /* void */
1072             }
1073         } /* 0 <= i < prMSeq->nseqs */
1074
1075     } /* there was a HMM batch file */
1076
1077
1078     if (0 < prOpts->iHMMInputFiles) {
1079         int iHMMInfileIndex;
1080
1081         /**
1082          * @warning old structure used to be initialised like this:
1083          * hmm_light rHMM = {0};
1084          */
1085         prHMMs = (hmm_light *) CKMALLOC( (prOpts->iHMMInputFiles) * sizeof(hmm_light));
1086
1087         for (iHMMInfileIndex=0; iHMMInfileIndex<prOpts->iHMMInputFiles; iHMMInfileIndex++) {
1088             char *pcHMMInput = prOpts->ppcHMMInput[iHMMInfileIndex];
1089             if (OK != readHMMWrapper(&prHMMs[iHMMInfileIndex], pcHMMInput)){
1090                 Log(&rLog, LOG_ERROR, "Processing of HMM file %s failed", pcHMMInput);
1091                 return -1;
1092             }
1093
1094 #if 0
1095             Log(&rLog, LOG_FORCED_DEBUG, "HMM length is %d", prHMMs[iHMMInfileIndex].L);
1096             Log(&rLog, LOG_FORCED_DEBUG, "n-display  is %d", prHMMs[iHMMInfileIndex].n_display);
1097             for (i = 0; NULL != prHMMs[prOpts->iHMMInputFiles].seq[i]; i++){
1098                 printf("seq[%d]: %s\n", i, prHMMs[iHMMInfileIndex].seq[i]);
1099             }
1100             Log(&rLog, LOG_FORCED_DEBUG, "Neff_HMM   is %f", prHMMs[iHMMInfileIndex].Neff_HMM);
1101 #endif
1102             if (rLog.iLogLevelEnabled <= LOG_DEBUG){
1103                 Log(&rLog, LOG_DEBUG, "print frequencies");
1104                 for (i = 0; i < prHMMs[iHMMInfileIndex].L; i++){
1105 #define PRINT_TAIL 5
1106                     if ( (PRINT_TAIL+1 == i) && (prHMMs[iHMMInfileIndex].L-PRINT_TAIL != i) ){
1107                         printf("....\n");
1108                     }
1109                     if ( (i > PRINT_TAIL) && (i < prHMMs[iHMMInfileIndex].L-PRINT_TAIL) ){
1110                         continue;
1111                     }
1112                     printf("%3d:", i);
1113                     for (j = 0; j < 20; j++){
1114                         printf("\t%1.3f", prHMMs[iHMMInfileIndex].f[i][j]);
1115                     }
1116                     printf("\n");
1117                 }
1118             } /* debug print block */
1119
1120             CKFREE(prOpts->ppcHMMInput[iHMMInfileIndex]);
1121         } /* for each background HMM file */
1122         CKFREE(prOpts->ppcHMMInput);
1123     } /* there were background HMM files */
1124
1125     /** read HMMs specific to individual sequences
1126      */
1127     if (iHMMbatch > 0){
1128         int i;
1129
1130         prHMMs = (hmm_light *) realloc( prHMMs, (prOpts->iHMMInputFiles + iHMMbatch + 1) * sizeof(hmm_light));
1131
1132         for (i = 0; i < iHMMbatch; i++){
1133             char *pcHMMInput = ppcHMMbatch[i];
1134
1135             if (OK != readHMMWrapper(&prHMMs[i + prOpts->iHMMInputFiles], pcHMMInput)){
1136                 Log(&rLog, LOG_ERROR, "Processing of HMM file %s failed", pcHMMInput);
1137                 return -1;
1138             }
1139
1140         } /* 0 <= i < iHMMbatch */
1141
1142     } /* there were HMM batch files */
1143
1144
1145     /* If the input ("non-profile") sequences are aligned, then turn
1146      * the alignment into a HMM and add to the list of background HMMs
1147      *
1148      */
1149     if (TRUE == prMSeq->aligned) {
1150         /* FIXME: gcc warns about missing initialiser here (-Wall -Wextra -pedantic) */
1151         hmm_light rHMMLocal = {0};
1152
1153         Log(&rLog, LOG_INFO,
1154             "Input sequences are aligned. Will turn alignment into HMM and add it to the user provided background HMMs.");
1155         /* certain gap parameters ('~' MSF) cause problems,
1156            sanitise them; FS, r258 -> r259 */
1157         SanitiseUnknown(prMSeq);
1158         if (OK !=
1159 #if INDIRECT_HMM
1160             AlnToHMM(&rHMMLocal, prMSeq)
1161 #else
1162             AlnToHMM2(&rHMMLocal, prOpts->rHhalignPara, prMSeq->seq, prMSeq->nseqs)
1163 #endif
1164             ) {
1165             Log(&rLog, LOG_ERROR, "Couldn't convert aligned input sequences to HMM. Will try to continue");
1166         } else {
1167             prHMMs = (hmm_light *) CKREALLOC(prHMMs, ((prOpts->iHMMInputFiles+1) * sizeof(hmm_light)));
1168             memcpy(&(prHMMs[prOpts->iHMMInputFiles]), &rHMMLocal, sizeof(hmm_light));
1169             prOpts->iHMMInputFiles++;
1170         }
1171     }
1172
1173
1174     /* If we have a profile turn it into a HMM and add to
1175      * the list of background HMMs.
1176      *
1177      */
1178     if (NULL != prMSeqProfile) {
1179         /* FIXME: gcc warns about missing initialiser here (-Wall -Wextra -pedantic) */
1180         hmm_light rHMMLocal = {0};
1181         Log(&rLog, LOG_INFO,
1182             "Turning profile1 into HMM and will use it during progressive alignment.");
1183         if (OK !=
1184 #if INDIRECT_HMM
1185             AlnToHMM(&rHMMLocal, prMSeqProfile)
1186 #else
1187             AlnToHMM2(&rHMMLocal, prOpts->rHhalignPara, prMSeqProfile->seq, prMSeqProfile->nseqs)
1188 #endif
1189             ) {
1190             Log(&rLog, LOG_ERROR, "Couldn't convert profile1 to HMM. Will try to continue");
1191         } else {
1192             prHMMs = (hmm_light *) CKREALLOC(prHMMs, ((prOpts->iHMMInputFiles+1) * sizeof(hmm_light)));
1193             memcpy(&(prHMMs[prOpts->iHMMInputFiles]), &rHMMLocal, sizeof(hmm_light));
1194             prOpts->iHMMInputFiles++;
1195         }
1196     }
1197
1198
1199     /* Now do a first alignment of the input sequences (prMSeq) adding
1200      * all collected background HMMs
1201      *
1202      */
1203     /* Determine progressive alignment order
1204      */
1205     if (TRUE == prMSeq->aligned) {
1206         if ( (SEQTYPE_PROTEIN == prMSeq->seqtype) && (TRUE == prOpts->bUseKimura) ){
1207             Log(&rLog, LOG_INFO, "%s %s",
1208                 "Input sequences are aligned.",
1209                 "Will use Kimura distances of aligned sequences.");
1210             prOpts->iPairDistType = PAIRDIST_SQUIDID_KIMURA;
1211         }
1212         else {
1213             prOpts->iPairDistType = PAIRDIST_SQUIDID;
1214         }
1215     }
1216
1217 #if 0
1218     Log(&rLog, LOG_WARN, "Using a sequential alignment order.");
1219     SequentialAlignmentOrder(&piOrderLR, prMSeq->nseqs);
1220 #else
1221     if (OK != AlignmentOrder(&piOrderLR, &pdSeqWeights, prMSeq,
1222                              prOpts->iPairDistType,
1223                              prOpts->pcDistmatInfile, prOpts->pcDistmatOutfile,
1224                              prOpts->iClusteringType, prOpts->iClustersizes,
1225                              prOpts->pcGuidetreeInfile, prOpts->pcGuidetreeOutfile, prOpts->pcClustfile,
1226                              prOpts->bUseMbed, prOpts->bPercID)) {
1227         Log(&rLog, LOG_ERROR, "AlignmentOrder() failed. Cannot continue");
1228         return -1;
1229     }
1230 #endif
1231
1232     /* if max-hmm-iter is set < 0 then do not perform alignment
1233      * there is a problem/feature(?) that the un-aligned sequences are output
1234      */
1235     if (prOpts->iMaxHMMIterations < 0){
1236         Log(&rLog, LOG_VERBOSE,
1237             "iMaxHMMIterations < 0 (%d), will not perform alignment", prOpts->iMaxHMMIterations);
1238         if (NULL != piOrderLR){
1239             free(piOrderLR); piOrderLR = NULL;
1240         }
1241         return 0;
1242     }
1243
1244
1245     /* Progressive alignment of input sequences. Order defined by
1246      * branching of guide tree (piOrderLR). Use optional
1247      * background HMM information (prHMMs[0..prOpts->iHMMInputFiles-1])
1248      *
1249      */
1250     dAlnScore = HHalignWrapper(prMSeq, piOrderLR, pdSeqWeights,
1251                                2*prMSeq->nseqs -1/* nodes */,
1252                                prHMMs, prOpts->iHMMInputFiles, -1, prOpts->rHhalignPara);
1253     dLastAlnScore = dAlnScore;
1254     Log(&rLog, LOG_VERBOSE,
1255         "Alignment score for first alignment = %f", dAlnScore);
1256
1257
1258
1259
1260     /* ------------------------------------------------------------
1261      *
1262      * prMSeq is aligned now. Now start iterations if requested and save the
1263      * alignment at the very end.
1264      *
1265      * @note We discard the background HMM information at this point,
1266      * because it was already used. Could consider to make this choice
1267      * optional.  FIXME
1268      *
1269      * ------------------------------------------------------------ */
1270
1271
1272     /* iteration after first alignment was computed (if not profile-profile
1273      * alignment)
1274      *
1275      */
1276     for (iIterationCounter=0;
1277          (iIterationCounter < prOpts->iNumIterations || prOpts->bIterationsAuto);
1278          iIterationCounter++) {
1279
1280         hmm_light rHMMLocal = {0};
1281         /* FIXME Keep copy of old alignment in case new one sucks? */
1282
1283
1284         if (iIterationCounter >= prOpts->iMaxHMMIterations
1285             &&
1286             iIterationCounter >= prOpts->iMaxGuidetreeIterations) {
1287             Log(&rLog, LOG_VERBOSE, "Reached maximum number of HMM and guide-tree iterations");
1288             break;
1289         }
1290
1291         if (! prOpts->bIterationsAuto) {
1292             Log(&rLog, LOG_INFO, "Iteration step %d out of %d",
1293                  iIterationCounter+1, prOpts->iNumIterations);
1294         } else {
1295             Log(&rLog, LOG_INFO, "Iteration step %d out of <auto>",
1296                  iIterationCounter+1);
1297         }
1298 #if 0
1299         if (rLog.iLogLevelEnabled <= LOG_VERBOSE) {
1300             char zcIntermediate[1000] = {0};
1301             char *pcFormat = "fasta";
1302             sprintf(zcIntermediate, "clustalo-aln-iter~%d~", iIterationCounter);
1303 #define LINE_WRAP 60
1304             if (WriteAlignment(prMSeq, zcIntermediate, MSAFILE_A2M, LINE_WRAP)) {
1305                 Log(&rLog, LOG_ERROR, "Could not save alignment to %s", zcIntermediate);
1306                 return -1;
1307             }
1308         }
1309 #endif
1310
1311
1312         /* new guide-tree
1313          *
1314          */
1315         if (iIterationCounter < prOpts->iMaxGuidetreeIterations) {
1316             /* determine progressive alignment order
1317              *
1318              * few things are different now when calling AlignmentOrder:
1319              * - we have to ignore prOpts->pcDistmatInfile and pcGuidetreeInfile
1320              *   as they were used before
1321              * - the corresponding outfiles are still valid though
1322              */
1323             /* Free stuff that has already been allocated by or further
1324              * downstream of AlignmentOrder()
1325              */
1326             if (NULL != piOrderLR)
1327                 CKFREE(piOrderLR);
1328             if (NULL != pdSeqWeights)
1329                 CKFREE(pdSeqWeights);
1330             Log(&rLog, LOG_INFO, "Computing new guide tree (iteration step %d)");
1331             if (AlignmentOrder(&piOrderLR, &pdSeqWeights, prMSeq,
1332                                ((SEQTYPE_PROTEIN == prMSeq->seqtype) && (TRUE == prOpts->bUseKimura)) ? PAIRDIST_SQUIDID_KIMURA : PAIRDIST_SQUIDID,
1333                                NULL, prOpts->pcDistmatOutfile,
1334                                prOpts->iClusteringType, prOpts->iClustersizes,
1335                                NULL, prOpts->pcGuidetreeOutfile, prOpts->pcClustfile,
1336                                prOpts->bUseMbedForIteration, prOpts->bPercID)) {
1337                 Log(&rLog, LOG_ERROR, "AlignmentOrder() failed. Cannot continue");
1338                 return -1;
1339             }
1340         } else {
1341             Log(&rLog, LOG_INFO, "Skipping guide-tree iteration at iteration step %d (reached maximum)",
1342                 iIterationCounter);
1343         }
1344
1345
1346         /* new local hmm iteration
1347          *
1348          */
1349         /* non-residue/gap characters will crash AlnToHMM2(),
1350            therefore sanitise unknown characters, FS, r259 -> r260 */
1351         SanitiseUnknown(prMSeq);
1352         if (iIterationCounter < prOpts->iMaxHMMIterations) {
1353             Log(&rLog, LOG_INFO, "Computing HMM from alignment");
1354
1355             if (OK !=
1356 #if INDIRECT_HMM
1357                 AlnToHMM(&rHMMLocal, prMSeq)
1358 #else
1359                 AlnToHMM2(&rHMMLocal, prOpts->rHhalignPara, prMSeq->seq, prMSeq->nseqs)
1360 #endif
1361                 ) {
1362                 Log(&rLog, LOG_ERROR, "Couldn't convert alignment to HMM. Will stop iterating now...");
1363                 break;
1364             }
1365         } else {
1366             Log(&rLog, LOG_INFO, "Skipping HMM iteration at iteration step %d (reached maximum)",
1367                  iIterationCounter);
1368         }
1369
1370
1371         /* align the sequences (again)
1372          */
1373         dAlnScore = HHalignWrapper(prMSeq, piOrderLR, pdSeqWeights,
1374                                    2*prMSeq->nseqs -1/* nodes */, &rHMMLocal, 1, -1,
1375                                    prOpts->rHhalignPara);
1376         Log(&rLog, LOG_VERBOSE,
1377             "Alignment score for alignmnent in hmm-iteration no %d = %f (last score = %f)",
1378              iIterationCounter+1, dAlnScore, dLastAlnScore);
1379
1380
1381         FreeHMMstruct(&rHMMLocal);
1382
1383 #if 0
1384         /* FIXME: need a better score for automatic iteration */
1385         if (prOpts->bIterationsAuto) {
1386             /* automatic iteration: break if score improvement was not
1387              * big enough
1388              */
1389             double dScoreImprovement = (dAlnScore-dLastAlnScore)/dLastAlnScore;
1390             if (dScoreImprovement < ITERATION_SCORE_IMPROVEMENT_THRESHOLD) {
1391                 Log(&rLog, LOG_INFO,
1392                     "Stopping after %d guide-tree iterations. No further alignment score improvement achieved.",
1393                      iIterationCounter+1);
1394                 /* use previous alignment */
1395                 FreeMSeq(&prMSeq);
1396                 Log(&rLog, LOG_FORCED_DEBUG, "FIXME: %s", "CopyMSeq breaks things in this context");
1397                 CopyMSeq(&prMSeq, prMSeqCopy);
1398                 /* FIXME: prOpts->pcDistmatOutfile and pcGuidetreeOutfile
1399                  * might have been updated, but then discarded here?
1400                  */
1401                 break;
1402             } else {
1403                 Log(&rLog, LOG_INFO,
1404                     "Got a %d%% better score in iteration step %d",
1405                      (int)dScoreImprovement*100, iIterationCounter+1);
1406                 FreeMSeq(&prMSeqCopy);
1407             }
1408         }
1409         dLastAlnScore = dAlnScore;
1410 #endif
1411
1412     }
1413     /* end of iterations */
1414
1415
1416
1417     /* Last step: if a profile was also provided then align now-aligned mseq
1418      * with this profile
1419      *
1420      * Don't use the backgrounds HMMs anymore and don't iterate.
1421      * (which was done before).
1422      *
1423      */
1424     if (NULL != prMSeqProfile) {
1425         if (AlignProfiles(prMSeq, prMSeqProfile, prOpts->rHhalignPara)) {
1426             Log(&rLog, LOG_ERROR, "An error occured during the profile/profile alignment");
1427             return -1;
1428         }
1429     }
1430
1431     if (NULL != prOpts->pcPosteriorfile){
1432
1433         hmm_light rHMMLocal = {0};
1434
1435         if (OK !=
1436 #if INDIRECT_HMM
1437             AlnToHMM(&rHMMLocal, prMSeq)
1438 #else
1439             AlnToHMM2(&rHMMLocal, prOpts->rHhalignPara, prMSeq->seq, prMSeq->nseqs)
1440 #endif
1441             ) {
1442             Log(&rLog, LOG_ERROR, "Couldn't convert alignment to HMM. Will not do posterior probabilities...");
1443         }
1444         PosteriorProbabilities(prMSeq, rHMMLocal, prOpts->rHhalignPara, prOpts->pcPosteriorfile);
1445         FreeHMMstruct(&rHMMLocal);
1446     }
1447
1448     if (NULL != piOrderLR) {
1449         CKFREE(piOrderLR);
1450     }
1451     if (NULL != pdSeqWeights) {
1452         CKFREE(pdSeqWeights);
1453     }
1454     if (0 < prOpts->iHMMInputFiles) {
1455         for (i=0; i<prOpts->iHMMInputFiles; i++) {
1456             FreeHMMstruct(&prHMMs[i]);
1457         }
1458         CKFREE(prHMMs);
1459     }
1460
1461     return 0;
1462 }
1463 /* end of Align() */
1464
1465
1466
1467
1468 /**
1469  * @brief Align two profiles, ie two sets of prealigned sequences. Already
1470  * aligned columns won't be changed.
1471  *
1472  * @param[out] prMSeqProfile1
1473  * First profile/aligned set of sequences. Merged alignment will be found in
1474  * here.
1475  * @param[in] prMSeqProfile2
1476  * First profile/aligned set of sequences
1477  * @param[in] rHhalignPara
1478  * FIXME
1479  *
1480  * @return 0 on success, -1 on failure
1481  *
1482  */
1483 int
1484 AlignProfiles(mseq_t *prMSeqProfile1,
1485               mseq_t *prMSeqProfile2, hhalign_para rHhalignPara) {
1486
1487     double dAlnScore;
1488
1489     /* number of seqs in first half of joined profile */
1490     int iProfProfSeparator = prMSeqProfile1->nseqs;
1491
1492     assert(TRUE == prMSeqProfile1->aligned);
1493     assert(TRUE == prMSeqProfile2->aligned);
1494
1495     Log(&rLog, LOG_INFO, "Performing profile/profile alignment");
1496
1497     /* Combine the available mseqs into prMSeq
1498      * which will be aligned afterwards.
1499      */
1500     JoinMSeqs(&prMSeqProfile1, prMSeqProfile2);
1501
1502
1503     /* set alignment flag explicitly to FALSE */
1504     prMSeqProfile1->aligned = FALSE;
1505
1506     dAlnScore = HHalignWrapper(prMSeqProfile1,
1507                                NULL, /* no order */
1508                                NULL, /* no weights */
1509                                3, /* nodes: root+2profiles */
1510                                NULL, 0 /* no bg-hmms */,
1511                                iProfProfSeparator, rHhalignPara);
1512
1513     Log(&rLog, LOG_VERBOSE, "Alignment score is = %f", dAlnScore);
1514
1515     return 0;
1516 }
1517 /* end of AlignProfiles() */
1518
1519 /**
1520  * @brief read pseudo-count 'fudge' parameters from file
1521  *
1522  * @param[out] rHhalignPara_p
1523  * structure that holds several hhalign parameters
1524  * @param[in] pcPseudoFile
1525  * name of file with pseudo-count information.
1526  * format must be collection of pairs of lines where one line specifies name of parameter
1527  * (gapb,gapd,gape,gapf,gapg,gaph,gapi,pca,pcb,pcc,gapbV,gapdV,gapeV,gapfV,gapgV,gaphV,gapiV,pcaV,pcbV,pccV)
1528  * followed by second line with the (double) value of this parameter.
1529  *
1530  * order of parameters is not fixed, not all parameters have to be set
1531  */
1532 int ReadPseudoCountParams(hhalign_para *rHhalignPara_p, char *pcPseudoFile){
1533
1534
1535     FILE *pfIn = NULL;
1536     char zcLine[1000] = {0};
1537     char zcFudge[1000] = {0};
1538     double dVal = 0.00;
1539     char *pcToken = NULL;
1540     char *pcEnd = NULL;
1541
1542     if (NULL == (pfIn = fopen(pcPseudoFile, "r"))){
1543         Log(&rLog, LOG_FATAL, "File %s with pseudo-count parameters does not exist", pcPseudoFile);
1544     }
1545     else {
1546       while (NULL != fgets(zcLine, 1000, pfIn)){
1547
1548             if (NULL == (pcToken = strtok(zcLine, " \040\t\n"))){
1549                 Log(&rLog, LOG_FATAL, "no token specifying pseudo-count parameter in file %s\n"
1550                        , pcPseudoFile
1551                        );
1552             }
1553             strcpy(zcFudge, pcToken);
1554
1555             if (NULL == fgets(zcLine, 1000, pfIn)){
1556             Log(&rLog, LOG_FATAL, "no line with parameter in file %s associated with %s\n"
1557                 , pcPseudoFile, zcFudge
1558                      );
1559             }
1560             else {
1561               if (NULL == (pcToken = strtok(zcLine, " \040\t\n"))){
1562               Log(&rLog, LOG_FATAL, "no token in 2nd line in file %s associated with %s\n"
1563                        , pcPseudoFile, zcFudge
1564                        );
1565               }
1566               else {
1567                 dVal = (double)strtod(pcToken, &pcEnd);
1568                 if ('\0' != *pcEnd){
1569                     Log(&rLog, LOG_FATAL, "token in file %s associated with %s not valid double (%s)\n"
1570                         , pcPseudoFile, zcFudge, pcToken
1571                         );
1572                 }
1573               }
1574             }
1575 #if 0
1576             printf("%s:%s:%d: read file %s: fudge = %s, val = %f\n"
1577                    , __FUNCTION__, __FILE__, __LINE__
1578                    , pcPseudoFile, zcFudge, dVal
1579                    );
1580 #endif
1581
1582             if      (0 == strcmp(zcFudge, "gapb")){
1583             rHhalignPara_p->gapb = dVal;
1584             }
1585             else if (0 == strcmp(zcFudge, "gapd")){
1586             rHhalignPara_p->gapd = dVal;
1587             }
1588             else if (0 == strcmp(zcFudge, "gape")){
1589             rHhalignPara_p->gape = dVal;
1590             }
1591             else if (0 == strcmp(zcFudge, "gapf")){
1592             rHhalignPara_p->gapf = dVal;
1593             }
1594             else if (0 == strcmp(zcFudge, "gapg")){
1595             rHhalignPara_p->gapg = dVal;
1596             }
1597             else if (0 == strcmp(zcFudge, "gaph")){
1598             rHhalignPara_p->gaph = dVal;
1599             }
1600             else if (0 == strcmp(zcFudge, "gapi")){
1601             rHhalignPara_p->gapi = dVal;
1602             }
1603             else if (0 == strcmp(zcFudge, "pca")){
1604             rHhalignPara_p->pca = dVal;
1605             }
1606             else if (0 == strcmp(zcFudge, "pcb")){
1607             rHhalignPara_p->pcb = dVal;
1608             }
1609             else if (0 == strcmp(zcFudge, "pcc")){
1610             rHhalignPara_p->pcc = dVal;
1611             }
1612             else if (0 == strcmp(zcFudge, "gapbV")){
1613             rHhalignPara_p->gapbV = dVal;
1614             }
1615             else if (0 == strcmp(zcFudge, "gapdV")){
1616             rHhalignPara_p->gapdV = dVal;
1617             }
1618             else if (0 == strcmp(zcFudge, "gapeV")){
1619             rHhalignPara_p->gapeV = dVal;
1620             }
1621             else if (0 == strcmp(zcFudge, "gapfV")){
1622             rHhalignPara_p->gapfV = dVal;
1623             }
1624             else if (0 == strcmp(zcFudge, "gapgV")){
1625             rHhalignPara_p->gapgV = dVal;
1626             }
1627             else if (0 == strcmp(zcFudge, "gaphV")){
1628             rHhalignPara_p->gaphV = dVal;
1629             }
1630             else if (0 == strcmp(zcFudge, "gapiV")){
1631             rHhalignPara_p->gapiV = dVal;
1632             }
1633             else if (0 == strcmp(zcFudge, "pcaV")){
1634             rHhalignPara_p->pcaV = dVal;
1635             }
1636             else if (0 == strcmp(zcFudge, "pcbV")){
1637             rHhalignPara_p->pcbV = dVal;
1638             }
1639             else if (0 == strcmp(zcFudge, "pccV")){
1640             rHhalignPara_p->pccV = dVal;
1641             }
1642             else {
1643             Log(&rLog, LOG_FATAL, "%s not a valid pseudo-count parameter\n"
1644                 "must be one of [gapb,gapd,gape,gapf,gapg,gaph,gapi,pca,pcb,pcc,gapbV,gapdV,gapeV,gapfV,gapgV,gaphV,gapiV,pcaV,pcbV,pccV]\n"
1645                 , zcFudge);
1646             } /* switched between parameters */
1647
1648       } /* while !EOF */
1649       fclose(pfIn); pfIn = NULL;
1650
1651     } /* there was a parameter file */
1652 #if 0
1653     printf("%s:%s:%d: gapb=%g,gapd=%g,gape=%g,gapf=%g,gapg=%g,gaph=%g,gapi=%g,pca=%g,pcb=%g,pcc=%g\n"
1654            , __FUNCTION__, __FILE__, __LINE__
1655            , rHhalignPara_p->gapb, rHhalignPara_p->gapd, rHhalignPara_p->gape, rHhalignPara_p->gapf, rHhalignPara_p->gapg, rHhalignPara_p->gaph, rHhalignPara_p->gapi, rHhalignPara_p->pca, rHhalignPara_p->pcb, rHhalignPara_p->pcc
1656            );
1657 #endif
1658 #if 0
1659     printf("%s:%s:%d: gapbV=%g,gapdV=%g,gapeV=%g,gapfV=%g,gapgV=%g,gaphV=%g,gapiV=%g,pcaV=%g,pcbV=%g,pccV=%g\n"
1660            , __FUNCTION__, __FILE__, __LINE__
1661            , rHhalignPara_p->gapbV, rHhalignPara_p->gapdV, rHhalignPara_p->gapeV, rHhalignPara_p->gapfV, rHhalignPara_p->gapgV, rHhalignPara_p->gaphV, rHhalignPara_p->gapiV, rHhalignPara_p->pcaV, rHhalignPara_p->pcbV, rHhalignPara_p->pccV
1662            );
1663 #endif
1664
1665
1666     return 0;
1667
1668 } /* this is the end of ReadPseudoCountParams() */
1669
1670
1671 /*
1672  * EOF clustal-omega.c
1673  */