comparison clustalomega/clustal-omega-0.2.0/src/clustal-omega.c @ 0:ff1768533a07

Migrated tool version 0.2 from old tool shed archive to new tool shed repository

0:ff1768533a07

/* -*- mode: c; tab-width: 4; c-basic-offset: 4;  indent-tabs-mode: nil -*- */

/*********************************************************************
 * Clustal Omega - Multiple sequence alignment
 *
 * Copyright (C) 2010 University College Dublin
 *
 * Clustal-Omega is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 *
 * This file is part of Clustal-Omega.
 *
 ********************************************************************/

/*
 *  RCS $Id: clustal-omega.c 241 2011-05-04 14:37:17Z fabian $
 */

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <assert.h>

#include "clustal-omega.h"
#include "hhalign/general.h"

/* The following comment block contains the frontpage/mainpage of the doxygen
 *  documentation. Please add some more info. FIXME add more
 */

/**
 *
 * @mainpage Clustal-Omega Documentation
 *
 * @section intro_sec Introduction
 *
 * For more information see http://www.clustal.org/
 *
 * @section api_section API
 *
 * @subsection example_prog_subsection An Example Program
 *
 * To use libclustalo you will have to include the clustal-omega.h header and
 * link against libclustalo. For linking against libclustalo you will have to
 * use a C++ compiler, no matter if your program was written in C or C++.
 *
 * First compile (no linking) your source (for an example see section "\ref
 * example_src_subsubsec"):
 *
 * @code
 * $ gcc -c -ansi -Wall clustalo-api-test.c
 * @endcode
 *
 * Then link against libclustalo (we recommend the use of pkg-config as
 * explained in \ref pkgconfig_subsubsec). Assuming Clustal Omega was installed
 * in system-wide default directory (e.g. /usr) just type:
 *
 * @code
 * $ g++  -ansi -Wall -o clustalo-api-test clustalo-api-test.o  -lclustalo
 * @endcode
 *
 * Voila! Now you have your own alignment program which can be run with
 *
 * @code
 * $ ./clustalo-api-test <your-sequence-input>
 * @endcode
 *  
 * It's best to use the same compiler that you used for compiling libclustal.
 * If libclustal was compiled with OpenMP support, you will have to use OpenMP
 * flags for you program as well.
 *
 *
 * @subsubsection pkgconfig_subsubsec Using pkg-config / Figuring out compiler flags
 *
 * Clustal Omega comes with support for <a
 * href="http://pkg-config.freedesktop.org">pkg-config</a>, which means you
 * can run
 *
 * @code
 * $ pkg-config --cflags --libs clustalo
 * @endcode
 *
 * to figure out cflags and library flags needed to compile and link against
 * libclustalo. This is especially handy if Clustal Omega was installed to a
 * non-standard directory.
 *  
 * You might have to change PKG_CONFIG_PATH. For example, if you used the prefix $HOME/local/ for
 * installation then you will first need to set PKG_CONFIG_PATH:
 *
 * @code
 * $ export PKG_CONFIG_PATH=$HOME/local/lib/pkgconfig
 * $ pkg-config --cflags --libs clustalo
 * @endcode
 *  
 *  
 * To compile your source use:
 *  
 * @code
 * $ export PKG_CONFIG_PATH=$HOME/local/lib/pkgconfig
 * $ gcc -c -ansi -Wall clustalo-api-test.c  $(pkg-config --cflags clustalo)
 * $ g++  -ansi -Wall -o clustalo-api-test clustalo-api-test.o $(pkg-config --libs clustalo)
 * @endcode
 *
 *
 * @subsubsection example_src_subsubsec Example Source Code
 *
 * @include "clustalo-api-test.c"
 *
 *
 */




/* FIXME: doc */
/* the following are temporary flags while the code is still under construction;
   had problems internalising hhmake, so as temporary crutch 
   write alignment to file and get external hmmer/hhmake via system call 
   to read alignment and convert into HMM
   All this will go, once hhmake is properly internalised */
#define INDIRECT_HMM 0 /* temp flag: (1) write aln to file, use system(hmmer/hhmake), (0) internal hhmake */
#define USEHMMER 1 /* temp flag: use system(hmmer) to build HMM */
#define USEHHMAKE (!USEHMMER) /* temp flag: use system(hhmake) to build HMM */


/* shuffle order of input sequences */
#define SHUFFLE_INPUT_SEQ_ORDER 0

/* sort input sequences by length */
#define SORT_INPUT_SEQS 0


int iNumberOfThreads;

/* broken, unused and lonely */
static const int ITERATION_SCORE_IMPROVEMENT_THRESHOLD = 0.01;


/**
 * @brief Print Long version information to pre-allocated char. 
 *
 * @note short version
 * information is equivalent to PACKAGE_VERSION
 *
 * @param[out] pcStr
 * char pointer to write to preallocated to hold iSize chars.
 * @param[in] iSize
 * size of pcStr
 */
void 
PrintLongVersion(char *pcStr, int iSize)
{
    snprintf(pcStr, iSize, "version %s; code-name '%s'; build date %s",
             PACKAGE_VERSION, PACKAGE_CODENAME, __DATE__);
}
/* end of PrintLongVersion() */



/**
 * @brief free aln opts members
 *
 */
void
FreeAlnOpts(opts_t *prAlnOpts) {
    if (NULL != prAlnOpts->pcGuidetreeInfile) {
        CKFREE(prAlnOpts->pcGuidetreeInfile);
    }
    if (NULL != prAlnOpts->pcGuidetreeOutfile) {
        CKFREE(prAlnOpts->pcGuidetreeOutfile);
    }
    if (NULL  != prAlnOpts->pcDistmatOutfile) {
        CKFREE(prAlnOpts->pcDistmatOutfile);
    }
    if (NULL != prAlnOpts->pcDistmatInfile) {
        CKFREE(prAlnOpts->pcDistmatInfile);
    }
}
/* end of FreeAlnOpts() */



/**
 * @brief Sets members of given user opts struct to default values
 *
 * @param[out] prOpts
 * User opt struct to initialise
 *
 */
void
SetDefaultAlnOpts(opts_t *prOpts) {
    prOpts->bAutoOptions = FALSE;
    
    prOpts->pcDistmatInfile = NULL;
    prOpts->pcDistmatOutfile = NULL;
	
    prOpts->iClusteringType = CLUSTERING_UPGMA;
    prOpts->iPairDistType = PAIRDIST_KTUPLE;
    prOpts->bUseMbed = FALSE;
    prOpts->bUseMbedForIteration = FALSE;
    prOpts->pcGuidetreeOutfile = NULL;
    prOpts->pcGuidetreeInfile = NULL;
	
    prOpts->ppcHMMInput = NULL;
    prOpts->iHMMInputFiles = 0;
		
    prOpts->iNumIterations = 0;
    prOpts->bIterationsAuto = FALSE;
    prOpts->iMaxGuidetreeIterations = INT_MAX;
    prOpts->iMaxHMMIterations = INT_MAX;
    prOpts->iMacRam = 2048; /* give 2GB to MAC algorithm. FS, r240 -> r241 */
 }
/* end of SetDefaultAlnOpts() */



/**
 * @brief Check logic of parsed user options. Will exit (call Log(&rLog, LOG_FATAL, ))
 * on Fatal logic error
 *
 * @param[in] prOpts
 * Already parsed user options
 * 
 */    
void
AlnOptsLogicCheck(opts_t *prOpts)
{
    /* guide-tree & distmat
     *
     */
    if (prOpts->pcDistmatInfile && prOpts->pcGuidetreeInfile) {
        Log(&rLog, LOG_FATAL, "Read distances *and* guide-tree from file doesn't make sense.");
    }

    if (prOpts->pcDistmatOutfile && prOpts->pcGuidetreeInfile) {
        Log(&rLog, LOG_FATAL, "Won't be able to save distances to file, because I got a guide-tree as input.");
    }

    /* combination of options that don't make sense when not iterating
     */
    if (prOpts->iNumIterations==0 && prOpts->bIterationsAuto != TRUE)  {
        
        if (prOpts->pcGuidetreeInfile && prOpts->pcGuidetreeOutfile) {
            Log(&rLog, LOG_FATAL, "Got a guide-tree as input and output which doesn't make sense when not iterating.");            
        }
        /*
          if (prOpts->pcGuidetreeInfile && prOpts->bUseMbed > 0) {
          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.");            
          }
        */
        if (prOpts->bUseMbedForIteration > 0) {
            Log(&rLog, LOG_FATAL, "No iteration requested, but mbed for iteration was set. Paranoia exit.");
        }        
    }
    
    if (prOpts->iMacRam < 512) {

        Log(&rLog, LOG_INFO, "Memory for MAC Algorithm quite low, Viterbi Algorithm may be triggered.");

        if (prOpts->iMacRam < 1) {

            Log(&rLog, LOG_WARN, "Viterbi Algorithm always turned on, increase MAC-RAM to turn on MAC.");
        }
    }

    return; 
}
/* end of AlnOptsLogicCheck() */


/**
 * @brief FIXME doc
 */
void
PrintAlnOpts(FILE *prFile, opts_t *prOpts)
{
    int iAux;


    /* keep in same order as struct */
    fprintf(prFile, "option: auto-options = %d\n", prOpts->bAutoOptions);
    fprintf(prFile, "option: distmat-infile = %s\n", 
            NULL != prOpts->pcDistmatInfile? prOpts->pcDistmatInfile: "(null)");
	fprintf(prFile, "option: distmat-outfile = %s\n",
            NULL != prOpts->pcDistmatOutfile? prOpts->pcDistmatOutfile: "(null)");
	fprintf(prFile, "option: clustering-type = %d\n", prOpts->iClusteringType);
	fprintf(prFile, "option: pair-dist-type = %d\n", prOpts->iPairDistType);
	fprintf(prFile, "option: use-mbed = %d\n", prOpts->bUseMbed);
	fprintf(prFile, "option: use-mbed-for-iteration = %d\n", prOpts->bUseMbedForIteration);
	fprintf(prFile, "option: guidetree-outfile = %s\n", 
            NULL != prOpts->pcGuidetreeOutfile? prOpts->pcGuidetreeOutfile: "(null)");
	fprintf(prFile, "option: guidetree-infile = %s\n",
            NULL != prOpts->pcGuidetreeInfile? prOpts->pcGuidetreeInfile: "(null)");
    for (iAux=0; iAux<prOpts->iHMMInputFiles; iAux++) {
        fprintf(prFile, "option: hmm-input no %d = %s\n", iAux, prOpts->ppcHMMInput[iAux]);
    }
	fprintf(prFile, "option: hmm-input-files = %d\n", prOpts->iHMMInputFiles);
	fprintf(prFile, "option: num-iterations = %d\n", prOpts->iNumIterations);
	fprintf(prFile, "option: iterations-auto = %d\n", prOpts->bIterationsAuto);
	fprintf(prFile, "option: max-hmm-iterations = %d\n", prOpts->iMaxHMMIterations);
	fprintf(prFile, "option: max-guidetree-iterations = %d\n", prOpts->iMaxGuidetreeIterations);
}
/* end of PrintAlnOpts() */



/**
 * @brief Returns major version of HMMER. Whichever hmmbuild version
 * is found first in your PATH will be used
 *
 * @return -1 on error, major hmmer version otherwise
 *
 */
int
HmmerVersion()
{
    char zcHmmerTestCall[] = "hmmbuild -h";
    FILE *fp = NULL;
    int iMajorVersion = 0;
    char zcLine[16384];

    if (NULL == (fp = popen(zcHmmerTestCall, "r"))) {
        Log(&rLog, LOG_ERROR, "Couldn't exec %s", zcHmmerTestCall);
        return -1;
    }
    while (fgets(zcLine, sizeof(zcLine), fp)) {
        char *pcLocate;
        if ((pcLocate = strstr(zcLine, "HMMER "))) {
            iMajorVersion = atoi(&pcLocate[6]);
            break;
        }
    }
    pclose(fp);

    return iMajorVersion;
}
/* end of HmmerVersion() */



/**
 * @brief Create a HHM file from aligned sequences
 *
 * @warning Should be eliminated in the future 
 * as building routine should not create intermediate files
 *
 * @param[in] prMSeq
 * Aligned mseq_t
 * @param[in] pcHMMOut
 * HMM output file name
 *
 * @return Non-zero on error
 *
 */
int
AlnToHHMFile(mseq_t *prMSeq, char *pcHMMOut)
{
    char *tmp_aln = NULL;
    int retcode = OK;

    assert(NULL!=prMSeq);
    assert(NULL!=pcHMMOut);

    if (FALSE == prMSeq->aligned) {
        Log(&rLog, LOG_ERROR, "Sequences need to be aligned to create an HMM");
        return FAILURE;
    }

    /* Convert alignment to a2m, and call hhmake
     *
     * can't be static templates, or mktemp fails (at least on os x
     * (with a bus error))
     * 
     * gcc says we should use mkstemp to avoid race conditions, 
     * but that returns a file descriptor, which is of no use to 
     * us  
     */
    /* NOTE: the following won't work on windows: missing /tmp/ */
    tmp_aln = CkStrdup("/tmp/clustalo_tmpaln_XXXXXX");
    if (NULL == mktemp(tmp_aln)) {
        Log(&rLog, LOG_ERROR, "Could not create temporary alignment filename");
        retcode = FAILURE;
        goto cleanup_and_return;
    }
    if (WriteAlignment(prMSeq, tmp_aln, MSAFILE_A2M)) {
        Log(&rLog, LOG_ERROR, "Could not save alignment to %s", tmp_aln);
        retcode = FAILURE;
        goto cleanup_and_return;
    }

    if (HHMake_Wrapper(tmp_aln, pcHMMOut)){
        Log(&rLog, LOG_ERROR, "Could not convert alignment %s into HHM", tmp_aln);
        retcode = FAILURE;
        goto cleanup_and_return;
    }


 cleanup_and_return:

    if (NULL != tmp_aln) {
        if (FileExists(tmp_aln)) {
            if (remove(tmp_aln)) {
                Log(&rLog, LOG_WARN, "Removing %s failed. Continuing anyway", tmp_aln);
            }
        }
        CKFREE(tmp_aln);
    }

    return retcode; 

} /* end of AlnToHHMFile() */



/**
 * @brief Create a HMM file from aligned sequences
 *
 * @warning Should be replaced in the future by some internal HMM
 * building routine that does not call external programs
 *
 * @param[in] prMSeq
 * Aligned mseq_t
 * @param[in] pcHMMOut
 * HMM output file name
 *
 * @return Non-zero on error
 *

 */
int
AlnToHMMFile(mseq_t *prMSeq, const char *pcHMMOut)
{
    char *tmp_aln = NULL;
    char *tmp_hmm = NULL; /* only needed for hmmer3 to hmmer2 conversion */
    char cmdbuf[16384];
    int iHmmerVersion = 0;
    int retcode = OK;
    
    assert(NULL!=prMSeq);
    assert(NULL!=pcHMMOut);

    if (FALSE == prMSeq->aligned) {
        Log(&rLog, LOG_ERROR, "Sequences need to be aligned to create an HMM");
        return FAILURE;
    }

    iHmmerVersion = HmmerVersion();
    if (2 != iHmmerVersion && 3 != iHmmerVersion) {
        Log(&rLog, LOG_ERROR, "Could not find suitable HMMER binaries");
        return FAILURE;
    }

    /* Convert alignment to stockholm, call hmmbuild and then
     * either hmmconvert (hmmer3) or hmmcalibrate (hmmer2)
     *
     * can't be static templates, or mktemp fails (at least on os x
     * (with a bus error))
     *
     * gcc says we should use mkstemp to avoid race conditions,
     * but that returns a file descriptor, which is of no use to
     * us
     */
    /* NOTE: the following won't work on windows: missing /tmp/ */
    tmp_aln = CkStrdup("/tmp/clustalo_tmpaln_XXXXXX");
    if (NULL == mktemp(tmp_aln)) {
        Log(&rLog, LOG_ERROR, "Could not create temporary alignment filename");
        retcode = FAILURE;
        goto cleanup_and_return;
    }
    if (WriteAlignment(prMSeq, tmp_aln, MSAFILE_STOCKHOLM)) {
        Log(&rLog, LOG_ERROR, "Could not save alignment to %s", tmp_aln);
        retcode = FAILURE;
        goto cleanup_and_return;
    }

    if (2 == iHmmerVersion) {
        sprintf(cmdbuf, "hmmbuild %s %s >/dev/null && hmmcalibrate %s >/dev/null",
                pcHMMOut, tmp_aln, pcHMMOut);
        if (system(cmdbuf)) {
            Log(&rLog, LOG_ERROR, "Command '%s' failed", cmdbuf);
            retcode = FAILURE;
            goto cleanup_and_return;
        }
    } else if (3 == iHmmerVersion) {
        /* NOTE: the following won't work on windows: missing /tmp/ */
        tmp_hmm = CkStrdup("/tmp/clustalo_tmphmm2_XXXXXX");
        if (NULL == mktemp(tmp_hmm)) {
            Log(&rLog, LOG_ERROR, "Could not create temporary hmm filename");
            retcode = FAILURE;
            goto cleanup_and_return;
        }
        sprintf(cmdbuf, "hmmbuild %s %s >/dev/null && hmmconvert -2 %s > %s",
                tmp_hmm, tmp_aln, tmp_hmm, pcHMMOut);
        if (system(cmdbuf)) {
            Log(&rLog, LOG_ERROR, "Command '%s' failed", cmdbuf);
            retcode = FAILURE;
            goto cleanup_and_return;
        }
    } else {
        CKFREE(tmp_aln);
        Log(&rLog, LOG_FATAL, "Internal error: Unknown Hmmer version %d", iHmmerVersion);
    }


 cleanup_and_return:

    if (NULL != tmp_aln) {
        if (FileExists(tmp_aln)) {
            if (remove(tmp_aln)) {
                Log(&rLog, LOG_WARN, "Removing %s failed. Continuing anyway", tmp_aln);
            }
        }
        CKFREE(tmp_aln);
    }
    if (NULL != tmp_hmm) {
        if (FileExists(tmp_hmm)) {
            if (remove(tmp_hmm)) {
                Log(&rLog, LOG_WARN, "Removing %s failed. Continuing anyway", tmp_hmm);
            }
        }
        CKFREE(tmp_hmm);
    }

    return retcode;
}
/* end of AlnToHMMFile() */



/**
 * @brief Convert a multiple sequence structure into a HMM
 *
 * @param[out] prHMM
 * Pointer to preallocted HMM which will be set here
 * @param[in] prMSeq
 * Pointer to an alignment
 *
 * @return 0 on error, non-0 otherwise
 *
 * @see AlnToHMMFile()
 * 
 */    
int
AlnToHMM(hmm_light *prHMM, mseq_t *prMSeq)
{
    char *pcHMM; /* temp hmm file */

    Log(&rLog, LOG_INFO,
        "Using HMMER version %d to calculate a new HMM.",
         HmmerVersion());
    /* FIXME replace all this with internal HMM computation (HHmake) */

    /**
     * @warning the following probably won't work on windows: missing
     * /tmp/. Should be ok on Cygwin though
     */
    pcHMM = CkStrdup("/tmp/clustalo-hmm-iter_XXXXXX");
    if (NULL == mktemp(pcHMM)) {
        Log(&rLog, LOG_ERROR, "Could not create temporary hmm filename");
        CKFREE(pcHMM);
        return FAILURE;
    }
    
    /* Create a HMM representing the current alignment
     */
#if USEHMMER
    if (AlnToHMMFile(prMSeq, pcHMM)) {
        Log(&rLog, LOG_ERROR, "AlnToHMMFile() on %s failed.", pcHMM);
        CKFREE(pcHMM);
        return FAILURE;
    }
#elif USEHHMAKE
    if (AlnToHHMFile(prMSeq, pcHMM)) {
        Log(&rLog, LOG_ERROR, "AlnToHHMFile() on %s failed.", pcHMM);
        CKFREE(pcHMM);
        return FAILURE;
    }
    /* Log(&rLog, LOG_FATAL, "Method to create HHM (HMM using hhmake) not installed yet"); */
#else
    Log(&rLog, LOG_FATAL, "Unknown method to create temporary HMM");
#endif
    
    /* Read HMM information
     */
    if (OK != readHMMWrapper(prHMM, pcHMM)){
        Log(&rLog, LOG_ERROR, "Processing of HMM file %s failed", pcHMM);
        CKFREE(pcHMM);
        return FAILURE;
    }

    if (remove(pcHMM)) {
        Log(&rLog, LOG_WARN, "Removing %s failed. Continuing anyway", pcHMM);
    }
    CKFREE(pcHMM);
        
    return OK; 
}
/* end of AlnToHMM() */



/**
 * @brief FIXME
 *
 */
void
InitClustalOmega(int iNumThreadsRequested)
{

#ifdef HAVE_OPENMP
    iNumberOfThreads = iNumThreadsRequested;
    omp_set_num_threads(iNumberOfThreads);
#else
    if (iNumThreadsRequested>1) {
        Log(&rLog, LOG_FATAL, "Cannot change number of threads to %d. %s was build without OpenMP support.",
              iNumThreadsRequested, PACKAGE_NAME);
    }
    iNumberOfThreads = 1; /* need to set this, even if build without support */
#endif

    Log(&rLog, LOG_INFO, "Using %d threads",
         iNumberOfThreads);

}
/* end of InitClustalOmega() */



/**
 * @brief Defines an alignment order, which adds sequences
 * sequentially, i.e. one at a time starting with seq 1 & 2
 *
 * @param[out] piOrderLR_p
 * order in which nodes/profiles are to be merged/aligned
 * @param[in] iNumSeq
 * Number of sequences
 *
 * @see TraverseTree()
 *
 */
void
SequentialAlignmentOrder(int **piOrderLR_p, int iNumSeq)
{
    unsigned int uNodes = iNumSeq*2-1;
    unsigned int uNodeCounter = 0;
    unsigned int uSeqCounter = 0;

    Log(&rLog, LOG_FATAL, "FIXME: Untested...");
    
    (*piOrderLR_p) = (int *)CKCALLOC(DIFF_NODE * uNodes, sizeof(int));
    /* loop over merge nodes, which have per definition even indices
     * and set up children which have odd indices
     */
    uSeqCounter = 0;
    for (uNodeCounter=iNumSeq; uNodeCounter<uNodes; uNodeCounter+=1) {
         unsigned int uLeftChildNodeIndex = uNodeCounter-1;
         unsigned int uRightChildNodeIndex = uNodeCounter-iNumSeq+1;
         unsigned int uParentNodeIndex = uNodeCounter+1;

         /* merge node setup */
        (*piOrderLR_p)[DIFF_NODE*uNodeCounter+LEFT_NODE] = uLeftChildNodeIndex;
        (*piOrderLR_p)[DIFF_NODE*uNodeCounter+RGHT_NODE] = uRightChildNodeIndex;
        (*piOrderLR_p)[DIFF_NODE*uNodeCounter+PRNT_NODE] = uParentNodeIndex;
        /* only setup left child if at first merge node, all other left childs
         * should be merge nodes that are already set up. also correct
         * left node number here.
         */
        if (uNodeCounter==iNumSeq) {
            (*piOrderLR_p)[DIFF_NODE*uNodeCounter+LEFT_NODE] = 0;

            (*piOrderLR_p)[0+LEFT_NODE] = 0;
            (*piOrderLR_p)[0+RGHT_NODE] = 0;
            (*piOrderLR_p)[0+PRNT_NODE] = uNodeCounter;
            uSeqCounter++;

            Log(&rLog, LOG_FORCED_DEBUG, "Set up first leaf with node counter %d: left=%d right=%d parent=%d",
                      0,
                      (*piOrderLR_p)[DIFF_NODE*uLeftChildNodeIndex+LEFT_NODE],
                      (*piOrderLR_p)[DIFF_NODE*uLeftChildNodeIndex+RGHT_NODE],
                      (*piOrderLR_p)[DIFF_NODE*uLeftChildNodeIndex+PRNT_NODE]);
        }
        Log(&rLog, LOG_FORCED_DEBUG, "Set up merge node with node counter %d: left=%d right=%d parent=%d",
                  uNodeCounter, (*piOrderLR_p)[DIFF_NODE*uNodeCounter+LEFT_NODE],
                  (*piOrderLR_p)[DIFF_NODE*uNodeCounter+RGHT_NODE],
                  (*piOrderLR_p)[DIFF_NODE*uNodeCounter+PRNT_NODE]);
        
        /* right child */
        (*piOrderLR_p)[DIFF_NODE*uRightChildNodeIndex+LEFT_NODE] = uSeqCounter;
        (*piOrderLR_p)[DIFF_NODE*uRightChildNodeIndex+RGHT_NODE] = uSeqCounter;
        (*piOrderLR_p)[DIFF_NODE*uRightChildNodeIndex+PRNT_NODE] = uNodeCounter;
        uSeqCounter++;

        Log(&rLog, LOG_FORCED_DEBUG, "Set up leaf with node counter %d: left=%d right=%d parent=%d",
                  uRightChildNodeIndex, (*piOrderLR_p)[DIFF_NODE*uRightChildNodeIndex+LEFT_NODE],
                  (*piOrderLR_p)[DIFF_NODE*uRightChildNodeIndex+RGHT_NODE],
                  (*piOrderLR_p)[DIFF_NODE*uRightChildNodeIndex+PRNT_NODE]);
    }
}
/* end of SequentialAlignmentOrder() */



/**
 * @brief Defines the alignment order by calculating a guide tree. In
 * a first-step pairwise distances will be calculated (or read from a
 * file). In a second step those distances will be clustered and a
 * guide-tree created. Steps 1 and 2 will be skipped if a guide-tree
 * file was given, in which case the guide-tree will be just read from
 * the file.
 *
 * @param[out] piOrderLR_p
 * order in which nodes/profiles are to be merged/aligned
 * @param[out] pdSeqWeights_p
 * Sequence weights
 * @param[out] pdSeqWeights_p
 * Sequence weights
 * @param[in] prMSeq
 * The sequences from which the alignment order is to be calculated
 * @param[in] iPairDistType
 * Method of pairwise distance comparison
 * @param[in] pcDistmatInfile
 * If not NULL distances will be read from this file instead of being
 * calculated
 * @param[in] pcDistmatOutfile
 * If not NULL computed pairwise distances will be written to this file
 * @param[in] iClusteringType
 * Clustering method to be used to cluster the pairwise distances
 * @param[in] pcGuidetreeInfile
 * If not NULL guidetree will be read from this file. Skips pairwise
 * distance and guidetree computation
 * @param[in] pcGuidetreeOutfile
 * If not NULL computed guidetree will be written to this file
 * @param[in] bUseMbed
 * If TRUE, fast mBed guidetree computation will be employed
 *
 * @return Non-zero on error
 *
 */
int
AlignmentOrder(int **piOrderLR_p, double **pdSeqWeights_p, mseq_t *prMSeq,
               int iPairDistType, char *pcDistmatInfile, char *pcDistmatOutfile,
               int iClusteringType, char *pcGuidetreeInfile, char *pcGuidetreeOutfile,
               bool bUseMbed)
{
    /* pairwise distance matrix (tmat in 1.83) */
    symmatrix_t *distmat = NULL;
    /* guide tree */
    tree_t *prTree = NULL;
    int i = 0;

    
    /* Shortcut for only two sequences: Do not compute k-tuple
     * distances. Use the same logic as in TraverseTree() to setup
     * piOrderLR_p. Changes there will have to be reflected here as
     * well. */
    if (2==prMSeq->nseqs) {
        Log(&rLog, LOG_VERBOSE,
            "Have only two sequences: No need to compute pairwise score and compute a tree.");
        
        (*piOrderLR_p) = (int*) CKMALLOC(DIFF_NODE * 3 * sizeof(int));
        (*piOrderLR_p)[DIFF_NODE*0+LEFT_NODE] = 0;
        (*piOrderLR_p)[DIFF_NODE*0+RGHT_NODE] = 0;
        (*piOrderLR_p)[DIFF_NODE*0+PRNT_NODE] = 0;

        (*piOrderLR_p)[DIFF_NODE*1+LEFT_NODE] = 1;
        (*piOrderLR_p)[DIFF_NODE*1+RGHT_NODE] = 1;
        (*piOrderLR_p)[DIFF_NODE*1+PRNT_NODE] = 1;

        /* root */
        (*piOrderLR_p)[DIFF_NODE*2+LEFT_NODE] = 0;
        (*piOrderLR_p)[DIFF_NODE*2+RGHT_NODE] = 1;
        (*piOrderLR_p)[DIFF_NODE*2+PRNT_NODE] = 2;

        /* Same logic as CalcClustalWeights(). Changes there will
           have to be reflected here as well. */
#if USE_WEIGHTS
         (*pdWeights_p) = (double *) CKMALLOC(uNodeCount * sizeof(double));
         (*pdWeights_p)[0] = 0.5;
         (*pdWeights_p)[1] = 0.5;
#endif
        
        return OK;
    }

        
    /* compute distance & guide tree, alternatively read distances or
     * guide tree from file
     *
     */
    if (NULL != pcGuidetreeInfile) {
        Log(&rLog, LOG_INFO, "Reading guide-tree from %s", pcGuidetreeInfile);
        if (GuideTreeFromFile(&prTree, prMSeq, pcGuidetreeInfile)) {
            Log(&rLog, LOG_ERROR, "Reading of guide tree %s failed.", pcGuidetreeInfile);
            return FAILURE;
        }

    } else {

        if (bUseMbed) {
            if (Mbed(&prTree, prMSeq, iPairDistType, pcGuidetreeOutfile)) {
                Log(&rLog, LOG_ERROR, "mbed execution failed.");
                return FAILURE;
            }
            Log(&rLog, LOG_INFO, "Guide-tree computation (mBed) done.");
            if (NULL != pcDistmatOutfile) {
                Log(&rLog, LOG_INFO,
                    "Ignoring request to write distance matrix (am in mBed mode)");
            }
        } else {

            if (PairDistances(&distmat, prMSeq, iPairDistType,
                              0, prMSeq->nseqs, 0, prMSeq->nseqs,
                              pcDistmatInfile, pcDistmatOutfile)) {
                Log(&rLog, LOG_ERROR, "Couldn't compute pair distances");
                return FAILURE;
            }

            /* clustering of distances to get guide tree
             */
            if (CLUSTERING_UPGMA == iClusteringType) {
                char **labels;
                labels = (char**) CKMALLOC(prMSeq->nseqs * sizeof(char*));
                for (i=0; i<prMSeq->nseqs; i++) {
                    labels[i] = prMSeq->sqinfo[i].name;
                }

                GuideTreeUpgma(&prTree, labels, distmat, pcGuidetreeOutfile);
                Log(&rLog, LOG_INFO, "Guide-tree computation done.");

                CKFREE(labels);
            } else {
                Log(&rLog, LOG_FATAL, "INTERNAL ERROR %s",
                      "clustering method should have been checked before");
            }
        }
    }

#if USE_WEIGHTS
    /* derive sequence weights from tree
     *
     */
    Log(&rLog, LOG_INFO, "Calculating sequence weights");
    CalcClustalWeights(pdSeqWeights_p, prTree);
    for (i = 0; i < GetLeafCount(prTree); i++) {
        Log(&rLog, LOG_VERBOSE,
            "Weight for seq no %d: %s = %f",
             i, prMSeq->sqinfo[i].name, (*pdSeqWeights_p)[i]);
    }
#else
    Log(&rLog, LOG_DEBUG, "Not using weights");
#endif

    
    /* define traversing order of tree
     *
     */
    TraverseTree(piOrderLR_p, prTree, prMSeq);
    if (rLog.iLogLevelEnabled <= LOG_DEBUG) {
        /* FIXME: debug only, FS */
        uint uNodeIndex;
        FILE *fp = LogGetFP(&rLog, LOG_INFO);
        Log(&rLog, LOG_DEBUG, "left/right order after tree traversal");
        for (uNodeIndex = 0; uNodeIndex < GetNodeCount(prTree); uNodeIndex++) {
            fprintf(fp, "%3d:\t%2d/%2d -> %d\n", i,
                   (*piOrderLR_p)[DIFF_NODE*uNodeIndex+LEFT_NODE],
                   (*piOrderLR_p)[DIFF_NODE*uNodeIndex+RGHT_NODE],
                   (*piOrderLR_p)[DIFF_NODE*uNodeIndex+PRNT_NODE]);
        }
    }

    FreeMuscleTree(prTree);
    FreeSymMatrix(&distmat);

#if 0
    Log(&rLog, LOG_FATAL, "DEBUG EXIT before leaving %s", __FUNCTION__);
#endif
    return OK;
}
/* end of AlignmentOrder() */



/**
 * @brief Set some options automatically based on number of sequences. Might
 * overwrite some user-set options.
 *
 * @param[out] prOpts
 * Pointer to alignment options structure
 * @param[in] iNumSeq
 * Number of sequences to align
 */
void
SetAutoOptions(opts_t *prOpts, int iNumSeq) {

    Log(&rLog, LOG_INFO,
        "Setting options automatically based on input sequence characteristics (might overwrite some of your options).");
    
    if (iNumSeq >= 10000) {
        if (0 != prOpts->iNumIterations) {
            Log(&rLog, LOG_INFO, "Auto settings: Disabling iterations.");
            prOpts->iNumIterations = 0;
        }
        if (FALSE == prOpts->bUseMbed) {
            Log(&rLog, LOG_INFO, "Auto settings: Enabling mBed.");
            prOpts->bUseMbed = TRUE;
        }
        
    } else if (iNumSeq >= 1000) {
        if (1 != prOpts->iNumIterations) {
            Log(&rLog, LOG_INFO, "Auto settings: Setting iteration to 1.");
            prOpts->iNumIterations = 1;
        }
        if (FALSE == prOpts->bUseMbed) {
            Log(&rLog, LOG_INFO, "Auto settings: Enabling mBed.");
            prOpts->bUseMbed = TRUE;
        }
        
        
    } else {
        if (0 != prOpts->iNumIterations) {
            Log(&rLog, LOG_INFO, "Auto settings: Disabling iterations.");
            prOpts->iNumIterations = 0;
        }
        if (TRUE == prOpts->bUseMbed) {
            Log(&rLog, LOG_INFO, "Auto settings: Disabling mBed.");
            prOpts->bUseMbed = FALSE;
        }
    }
}
/* end of */



/**
 * @brief The main alignment function which wraps everything else.
 *
 * @param[out] prMSeq
 * *the* multiple sequences structure
 * @param[in] prMSeqProfile
 * optional profile to align against
 * @param[in] prOpts
 * alignmemnt options to use
 *
 * @return 0 on success, -1 on failure
 *
 */
int
Align(mseq_t *prMSeq, 
      mseq_t *prMSeqProfile,
      opts_t *prOpts,
      hhalign_para rHhalignPara) {
   
    /* HMM
     */
    /* structs with pseudocounts etc; one for each HMM infile, i.e.
     * index range: 0..iHMMInputFiles */
    hmm_light *prHMMs = NULL;

    /* MSA order in which nodes/profiles are to be merged/aligned
       (order of nodes in guide tree (left/right)*/
    int *piOrderLR = NULL;

    /* weights per sequence */
    double *pdSeqWeights = NULL;

    /* Iteration
     */
    int iIterationCounter = 0;
    double dAlnScore;
    /* last dAlnScore for iteration */
    double dLastAlnScore = -666.666;

    int i, j; /* aux */

    assert(NULL != prMSeq);
    if (NULL != prMSeqProfile) {
        assert(TRUE == prMSeqProfile->aligned);
    }


    /* automatic setting of options
     *
     */
    if (prOpts->bAutoOptions) {
        SetAutoOptions(prOpts, prMSeq->nseqs);
    }


#if SHUFFLE_INPUT_SEQ_ORDER
    /* 
     * shuffle input:  only useful for testing/debugging 
     */
    Log(&rLog, LOG_WARN, "Shuffling input sequences! (Will also change output order)");
    ShuffleMSeq(prMSeq);
#endif
        

#if SORT_INPUT_SEQS
    /* 
     * sort input:
     *
     * would ensure we *always* (unless we get into the mbed k-means stage)
     * get the same answer. usually you don't, because most pairwise alignment
     * scores are in theory not symmetric, therefore sequence ordering might
     * have an effect on the guide-tree. Sorting by length should get rid of
     * this (and takes no time even for 100k seqs). Benchmark results on
     * Balibase show almost no difference after sorting.
     */
    Log(&rLog, LOG_WARN, "Sorting input seq by length! This will also change the output order");
    SortMSeqByLength(prMSeq, 'd');

#endif


    /* Read backgrounds HMMs and store in prHMMs
     *
     */
    if (0 < prOpts->iHMMInputFiles) {
        int iHMMInfileIndex;
        
        /**
         * @warning old structure used to be initialised like this:
         * hmm_light rHMM = {0};
         */
        prHMMs = (hmm_light *) CKMALLOC(prOpts->iHMMInputFiles * sizeof(hmm_light));
        
        for (iHMMInfileIndex=0; iHMMInfileIndex<prOpts->iHMMInputFiles; iHMMInfileIndex++) {
            char *pcHMMInput = prOpts->ppcHMMInput[iHMMInfileIndex];
            if (OK != readHMMWrapper(&prHMMs[iHMMInfileIndex], pcHMMInput)){
                Log(&rLog, LOG_ERROR, "Processing of HMM file %s failed", pcHMMInput);
                return -1;
            }
            
#if 0
            Log(&rLog, LOG_FORCED_DEBUG, "HMM length is %d", prHMMs[iHMMInfileIndex].L);
            Log(&rLog, LOG_FORCED_DEBUG, "n-display  is %d", prHMMs[iHMMInfileIndex].n_display);
            for (i = 0; NULL != prHMMs[prOpts->iHMMInputFiles].seq[i]; i++){
                printf("seq[%d]: %s\n", i, prHMMs[iHMMInfileIndex].seq[i]);
            }
            Log(&rLog, LOG_FORCED_DEBUG, "Neff_HMM   is %f", prHMMs[iHMMInfileIndex].Neff_HMM);
#endif
            if (rLog.iLogLevelEnabled <= LOG_DEBUG){
                Log(&rLog, LOG_DEBUG, "print frequencies");
                for (i = 0; i < prHMMs[iHMMInfileIndex].L; i++){
#define PRINT_TAIL 5
                    if ( (PRINT_TAIL+1 == i) && (prHMMs[iHMMInfileIndex].L-PRINT_TAIL != i) ){
                        printf("....\n");
                    }
                    if ( (i > PRINT_TAIL) && (i < prHMMs[iHMMInfileIndex].L-PRINT_TAIL) ){
                        continue;
                    }
                    printf("%3d:", i);
                    for (j = 0; j < 20; j++){
                        printf("\t%1.3f", prHMMs[iHMMInfileIndex].f[i][j]);
                    }
                    printf("\n");
                }
            } /* debug print block */
            
            CKFREE(prOpts->ppcHMMInput[iHMMInfileIndex]);
        } /* for each background HMM file */
        CKFREE(prOpts->ppcHMMInput);
    } /* there were background HMM files */
    

    
    /* If the input ("non-profile") sequences are aligned, then turn
     * the alignment into a HMM and add to the list of background HMMs
     *
     */
    if (TRUE == prMSeq->aligned) {
        /* FIXME: gcc warns about missing initialiser here (-Wall -Wextra -pedantic) */
        hmm_light rHMMLocal = {0};
        
        Log(&rLog, LOG_INFO,
            "Input sequences are aligned. Will turn alignment into HMM and add it to the user provided background HMMs.");
        if (OK !=
#if INDIRECT_HMM 
            AlnToHMM(&rHMMLocal, prMSeq)
#else
            AlnToHMM2(&rHMMLocal, prMSeq->seq, prMSeq->nseqs)
#endif
            ) {
            Log(&rLog, LOG_ERROR, "Couldn't convert aligned input sequences to HMM. Will try to continue");
        } else {
            prHMMs = (hmm_light *) CKREALLOC(prHMMs, ((prOpts->iHMMInputFiles+1) * sizeof(hmm_light)));
            memcpy(&(prHMMs[prOpts->iHMMInputFiles]), &rHMMLocal, sizeof(hmm_light));
            prOpts->iHMMInputFiles++;
        }
    }

    
    /* If we have a profile turn it into a HMM and add to
     * the list of background HMMs.
     *
     */
    if (NULL != prMSeqProfile) {
        /* FIXME: gcc warns about missing initialiser here (-Wall -Wextra -pedantic) */
        hmm_light rHMMLocal = {0};
        Log(&rLog, LOG_INFO,
            "Turning profile1 into HMM and will use it during progressive alignment.");
        if (OK !=
#if INDIRECT_HMM 
            AlnToHMM(&rHMMLocal, prMSeqProfile)
#else 
            AlnToHMM2(&rHMMLocal, prMSeqProfile->seq, prMSeqProfile->nseqs)
#endif
            ) {
            Log(&rLog, LOG_ERROR, "Couldn't convert profile1 to HMM. Will try to continue");
        } else {
            prHMMs = (hmm_light *) CKREALLOC(prHMMs, ((prOpts->iHMMInputFiles+1) * sizeof(hmm_light)));
            memcpy(&(prHMMs[prOpts->iHMMInputFiles]), &rHMMLocal, sizeof(hmm_light));
            prOpts->iHMMInputFiles++;
        }
    }


    /* Now do a first alignment of the input sequences (prMSeq) adding
     * all collected background HMMs
     *
     */
    /* Determine progressive alignment order
     */
    if (TRUE == prMSeq->aligned) {
        Log(&rLog, LOG_INFO, "%s %s",
             "Input sequences are aligned.",
             "Will use Kimura distances of aligned sequences.");
        prOpts->iPairDistType = PAIRDIST_SQUIDID_KIMURA;
    }

#if 0
    Log(&rLog, LOG_WARN, "Using a sequential alignment order.");
    SequentialAlignmentOrder(&piOrderLR, prMSeq->nseqs);
#else
    if (OK != AlignmentOrder(&piOrderLR, &pdSeqWeights, prMSeq,
                             prOpts->iPairDistType,
                             prOpts->pcDistmatInfile, prOpts->pcDistmatOutfile,
                             prOpts->iClusteringType,
                             prOpts->pcGuidetreeInfile, prOpts->pcGuidetreeOutfile,
                             prOpts->bUseMbed)) {
        Log(&rLog, LOG_ERROR, "AlignmentOrder() failed. Cannot continue");
        return -1;
    }
#endif

    /* Progressive alignment of input sequences. Order defined by
     * branching of guide tree (piOrderLR). Use optional
     * background HMM information (prHMMs[0..prOpts->iHMMInputFiles-1])
     *
     */
    dAlnScore = HHalignWrapper(prMSeq, piOrderLR, pdSeqWeights,
                               2*prMSeq->nseqs -1/* nodes */,
                               prHMMs, prOpts->iHMMInputFiles, -1, rHhalignPara);
    dLastAlnScore = dAlnScore;
    Log(&rLog, LOG_VERBOSE,
        "Alignment score for first alignment = %f", dAlnScore);        




    /* ------------------------------------------------------------
     *
     * prMSeq is aligned now. Now start iterations if requested and save the
     * alignment at the very end.
     *
     * @note We discard the background HMM information at this point,
     * because it was already used. Could consider to make this choice
     * optional.  FIXME
     *
     * ------------------------------------------------------------ */

    
    /* iteration after first alignment was computed (if not profile-profile
     * alignment)
     *
     */
    for (iIterationCounter=0;
         (iIterationCounter < prOpts->iNumIterations || prOpts->bIterationsAuto);
         iIterationCounter++) {

        hmm_light rHMMLocal = {0};
        /* FIXME Keep copy of old alignment in case new one sucks? */


        if (iIterationCounter >= prOpts->iMaxHMMIterations
            &&
            iIterationCounter >= prOpts->iMaxGuidetreeIterations) {
            Log(&rLog, LOG_VERBOSE, "Reached maximum number of HMM and guide-tree iterations");
            break;
        }
        
        if (! prOpts->bIterationsAuto) {
            Log(&rLog, LOG_INFO, "Iteration step %d out of %d", 
                 iIterationCounter+1, prOpts->iNumIterations);
        } else {
            Log(&rLog, LOG_INFO, "Iteration step %d out of <auto>", 
                 iIterationCounter+1);
        }
#if 0
        if (rLog.iLogLevelEnabled <= LOG_VERBOSE) {
            char zcIntermediate[1000] = {0};
            char *pcFormat = "fasta";
            sprintf(zcIntermediate, "clustalo-aln-iter~%d~", iIterationCounter);
            if (WriteAlignment(prMSeq, zcIntermediate, MSAFILE_A2M)) {
                Log(&rLog, LOG_ERROR, "Could not save alignment to %s", zcIntermediate);
                return -1;
            }
        }
#endif


        /* new guide-tree
         *
         */
        if (iIterationCounter < prOpts->iMaxGuidetreeIterations) {
            /* determine progressive alignment order
             *
             * few things are different now when calling AlignmentOrder:
             * - we have to ignore prOpts->pcDistmatInfile and pcGuidetreeInfile
             *   as they were used before
             * - the corresponding outfiles are still valid though
             */
            /* Free stuff that has already been allocated by or further
             * downstream of AlignmentOrder()
             */
            if (NULL != piOrderLR)
                CKFREE(piOrderLR);
            if (NULL != pdSeqWeights)
                CKFREE(pdSeqWeights);
            if (AlignmentOrder(&piOrderLR, &pdSeqWeights, prMSeq,
                               PAIRDIST_SQUIDID_KIMURA /* override */, NULL, prOpts->pcDistmatOutfile,
                               prOpts->iClusteringType, NULL, prOpts->pcGuidetreeOutfile,
                               prOpts->bUseMbedForIteration)) {
                Log(&rLog, LOG_ERROR, "AlignmentOrder() failed. Cannot continue");
                return -1;
            }
        } else {
            Log(&rLog, LOG_INFO, "Skipping guide-tree iteration at iteration step %d (reached maximum)", 
                iIterationCounter);
        }


        /* new local hmm iteration
         *
         */
        if (iIterationCounter < prOpts->iMaxHMMIterations) {
            if (OK !=
#if INDIRECT_HMM 
                AlnToHMM(&rHMMLocal, prMSeq)
#else
                AlnToHMM2(&rHMMLocal, prMSeq->seq, prMSeq->nseqs)
#endif
                ) {
                Log(&rLog, LOG_ERROR, "Couldn't convert alignment to HMM. Will stop iterating now...");
                break;
            }
        } else {
            Log(&rLog, LOG_INFO, "Skipping HMM iteration at iteration step %d (reached maximum)", 
                 iIterationCounter);
        }

        
        /* align the sequences (again)
         */
        dAlnScore = HHalignWrapper(prMSeq, piOrderLR, pdSeqWeights,
                                   2*prMSeq->nseqs -1/* nodes */, &rHMMLocal, 1, -1, rHhalignPara);
        Log(&rLog, LOG_VERBOSE,
            "Alignment score for alignmnent in hmm-iteration no %d = %f (last score = %f)",
             iIterationCounter+1, dAlnScore, dLastAlnScore);
        
        
        FreeHMMstruct(&rHMMLocal);

#if 0
        /* FIXME: need a better score for automatic iteration */
        if (prOpts->bIterationsAuto) {
            /* automatic iteration: break if score improvement was not
             * big enough
             */
            double dScoreImprovement = (dAlnScore-dLastAlnScore)/dLastAlnScore;
            if (dScoreImprovement < ITERATION_SCORE_IMPROVEMENT_THRESHOLD) {
                Log(&rLog, LOG_INFO,
                    "Stopping after %d guide-tree iterations. No further alignment score improvement achieved.",
                     iIterationCounter+1);
                /* use previous alignment */
                FreeMSeq(&prMSeq);
                Log(&rLog, LOG_FORCED_DEBUG, "FIXME: %s", "CopyMSeq breaks things in this context");
                CopyMSeq(&prMSeq, prMSeqCopy);
                /* FIXME: prOpts->pcDistmatOutfile and pcGuidetreeOutfile
                 * might have been updated, but then discarded here?
                 */
                break;
            } else {
                Log(&rLog, LOG_INFO,
                    "Got a %d%% better score in iteration step %d",
                     (int)dScoreImprovement*100, iIterationCounter+1);
                FreeMSeq(&prMSeqCopy);
            }
        }
        dLastAlnScore = dAlnScore;
#endif

    }
    /* end of iterations */



    /* Last step: if a profile was also provided then align now-aligned mseq
     * with this profile
     *
     * Don't use the backgrounds HMMs anymore and don't iterate.
     * (which was done before).
     *
     */
    if (NULL != prMSeqProfile) {
        if (AlignProfiles(prMSeq, prMSeqProfile, rHhalignPara)) {
            Log(&rLog, LOG_ERROR, "An error occured during the profile/profile alignment");
            return -1;
        }
    }

    
    if (NULL != piOrderLR) {
        CKFREE(piOrderLR);
    }
    if (NULL != pdSeqWeights) {
        CKFREE(pdSeqWeights);
    }
    if (0 < prOpts->iHMMInputFiles) {
        for (i=0; i<prOpts->iHMMInputFiles; i++) {
            FreeHMMstruct(&prHMMs[i]);
        }
        CKFREE(prHMMs);
    }

    return 0;
}
/* end of Align() */




/**
 * @brief Align two profiles, ie two sets of prealigned sequences. Already
 * aligned columns won't be changed.
 *
 * @param[out] prMSeqProfile1
 * First profile/aligned set of sequences. Merged alignment will be found in
 * here.
 * @param[in] prMSeqProfile2
 * First profile/aligned set of sequences
 *
 * @return 0 on success, -1 on failure
 *
 */
int
AlignProfiles(mseq_t *prMSeqProfile1, 
      mseq_t *prMSeqProfile2, hhalign_para rHhalignPara) {
   
    double dAlnScore;

    /* number of seqs in first half of joined profile */
    int iProfProfSeparator = prMSeqProfile1->nseqs;

    assert(TRUE == prMSeqProfile1->aligned);
    assert(TRUE == prMSeqProfile2->aligned);

    Log(&rLog, LOG_INFO, "Performing profile/profile alignment");

    /* Combine the available mseqs into prMSeq
     * which will be aligned afterwards.
     */
    JoinMSeqs(&prMSeqProfile1, prMSeqProfile2);

        
    /* set alignment flag explicitly to FALSE */
    prMSeqProfile1->aligned = FALSE;
                    
    dAlnScore = HHalignWrapper(prMSeqProfile1,
                               NULL, /* no order */
                               NULL, /* no weights */
                               3, /* nodes: root+2profiles */
                               NULL, 0 /* no bg-hmms */,
                               iProfProfSeparator, rHhalignPara);
        
    Log(&rLog, LOG_VERBOSE, "Alignment score is = %f", dAlnScore);

    return 0;
}
/* end of AlignProfiles() */