new file mode 100644

@@ -0,0 +1,1435 @@

+/* -*- 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 284 2013-06-12 10:10:11Z fabian $

+ */

+

+#ifdef HAVE_CONFIG_H

+#include "config.h"

+#endif

+

+#include <assert.h>

+

+#include "clustal-omega.h"

+#include "hhalign/general.h"

+#include "clustal/hhalign_wrapper.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++. See

+ * below (\ref pkgconfig_subsubsec)) on how to figure out compiler flags with

+ * pkg-config.

+ *

+ * Assuming Clustal Omega was installed in system-wide default directory (e.g.

+ * /usr), first compile (don't link yet) your source (for example code see

+ * section \ref example_src_subsubsec) and then link against libclustalo:

+ *

+ * @code

+ * $ gcc -c -ansi -Wall clustalo-api-test.c

+ * $ g++  -ansi -Wall -o clustalo-api-test clustalo-api-test.o  -lclustalo

+ * @endcode

+ *

+ * Voila! Now you have your own alignment program based on Clustal Omega 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 as above but this time using proper flags:

+ *  

+ * @code

+ * $ export PKG_CONFIG_PATH=$HOME/local/lib/pkgconfig

+ * $ gcc -c -ansi -Wall $(pkg-config --cflags clustalo) clustalo-api-test.c  

+ * $ g++  -ansi -Wall -o clustalo-api-test $(pkg-config --libs clustalo) clustalo-api-test.o 

+ * @endcode

+ *

+ *

+ * @subsubsection example_src_subsubsec Example Source Code

+ *

+ * @include "clustalo-api-test.c"

+ *

+ *

+ */

+

+

+/* 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->iClustersizes = 100; /* FS, r274 -> */

+    prOpts->pcClustfile = NULL; /* FS, r274 -> */

+    prOpts->iClusteringType = CLUSTERING_UPGMA;

+    prOpts->iPairDistType = PAIRDIST_KTUPLE;

+    prOpts->bUseMbed = TRUE; /* FS, r250 -> */

+    prOpts->bUseMbedForIteration = TRUE; /* FS, r250 -> */

+    prOpts->pcGuidetreeOutfile = NULL;

+    prOpts->pcGuidetreeInfile = NULL;

+    prOpts->bPercID = FALSE;

+

+    prOpts->ppcHMMInput = NULL;

+    prOpts->iHMMInputFiles = 0;

+		

+    prOpts->iNumIterations = 0;

+    prOpts->bIterationsAuto = FALSE;

+    prOpts->iMaxGuidetreeIterations = INT_MAX;

+    prOpts->iMaxHMMIterations = INT_MAX;

+

+    SetDefaultHhalignPara(& prOpts->rHhalignPara);

+ }

+/* 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.");            

+          }

+        */

+        /*

+          AW: bUseMbedForIteration default since at least R252

+          if (prOpts->bUseMbedForIteration > 0) {

+            Log(&rLog, LOG_FATAL, "No iteration requested, but mbed for iteration was set. Paranoia exit.");

+          }        

+        */

+    }

+    

+    if (prOpts->rHhalignPara.iMacRamMB < 512) {

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

+        if (prOpts->rHhalignPara.iMacRamMB < 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);

+    fprintf(prFile, "option: iMacRamMB = %d\n", prOpts->rHhalignPara.iMacRamMB);

+    fprintf(prFile, "option: percent-id = %d\n", prOpts->bPercID);

+    fprintf(prFile, "option: use-kimura = %d\n", prOpts->bUseKimura);

+

+}

+/* 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;

+

+    int msa_c = 0;

+    char **msa_v = NULL;

+

+    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;

+    }

+#define LINE_WRAP 60

+

+    if (WriteAlignment(prMSeq, tmp_aln, MSAFILE_A2M, LINE_WRAP, FALSE, &msa_c, &msa_v)) {

+        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;

+    

+    int msa_c = 0;

+    char **msa_v = NULL;

+

+    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;

+    }

+#define LINE_WRAP 60

+    if (WriteAlignment(prMSeq, tmp_aln, MSAFILE_STOCKHOLM, LINE_WRAP, FALSE, &msa_c, &msa_v)) {

+        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 WIN32

+    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 */

+#elif 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, int iClustersizes, 

+               char *pcGuidetreeInfile, char *pcGuidetreeOutfile, char *pcClusterFile, 

+               bool bUseMbed, bool bPercID)

+{

+    /* 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.");

+

+        if (NULL != pcDistmatOutfile){

+            Log(&rLog, LOG_WARN, "Have only two sequences: Will not calculate/print distance matrix.");

+        }

+

+        (*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, iClustersizes, pcClusterFile)) {

+                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, bPercID, 

+                              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");

+            }

+        } /* did not use mBed */

+    } /* had to calculate tree (did not read from file) */

+

+

+#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