| dafeef0b |
/* -*- 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
*
*/
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).");
/* AW: new version of mbed is always good (uses subclusters) */
if (FALSE == prOpts->bUseMbed) {
Log(&rLog, LOG_INFO, "Auto settings: Enabling mBed.");
prOpts->bUseMbed = TRUE;
}
if (iNumSeq >= 1000) {
if (0 != prOpts->iNumIterations) {
Log(&rLog, LOG_INFO, "Auto settings: Disabling iterations.");
prOpts->iNumIterations = 0;
}
} else if (iNumSeq < 1000) {
if (1 != prOpts->iNumIterations) {
Log(&rLog, LOG_INFO, "Auto settings: Setting iteration to 1.");
prOpts->iNumIterations = 1;
}
}
}
/* 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
* alignment options to use
*
* @return 0 on success, -1 on failure
*
*/
int
Align(mseq_t *prMSeq,
mseq_t *prMSeqProfile,
opts_t *prOpts) {
/* 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.");
/* certain gap parameters ('~' MSF) cause problems,
sanitise them; FS, r258 -> r259 */
SanitiseUnknown(prMSeq);
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) {
if ( (SEQTYPE_PROTEIN == prMSeq->seqtype) && (TRUE == prOpts->bUseKimura) ){
Log(&rLog, LOG_INFO, "%s %s",
"Input sequences are aligned.",
"Will use Kimura distances of aligned sequences.");
prOpts->iPairDistType = PAIRDIST_SQUIDID_KIMURA;
}
else {
prOpts->iPairDistType = PAIRDIST_SQUIDID;
}
}
#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->iClustersizes,
prOpts->pcGuidetreeInfile, prOpts->pcGuidetreeOutfile, prOpts->pcClustfile,
prOpts->bUseMbed, prOpts->bPercID)) {
Log(&rLog, LOG_ERROR, "AlignmentOrder() failed. Cannot continue");
return -1;
}
#endif
/* if max-hmm-iter is set < 0 then do not perform alignment
* there is a problem/feature(?) that the un-aligned sequences are output
*/
if (prOpts->iMaxHMMIterations < 0){
Log(&rLog, LOG_VERBOSE,
"iMaxHMMIterations < 0 (%d), will not perform alignment", prOpts->iMaxHMMIterations);
return 0;
}
/* 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, prOpts->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);
#define LINE_WRAP 60
if (WriteAlignment(prMSeq, zcIntermediate, MSAFILE_A2M, LINE_WRAP, &msa_c, &msa_v)) {
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);
Log(&rLog, LOG_INFO, "Computing new guide tree (iteration step %d)");
if (AlignmentOrder(&piOrderLR, &pdSeqWeights, prMSeq,
((SEQTYPE_PROTEIN == prMSeq->seqtype) && (TRUE == prOpts->bUseKimura)) ? PAIRDIST_SQUIDID_KIMURA : PAIRDIST_SQUIDID,
NULL, prOpts->pcDistmatOutfile,
prOpts->iClusteringType, prOpts->iClustersizes,
NULL, prOpts->pcGuidetreeOutfile, prOpts->pcClustfile,
prOpts->bUseMbedForIteration, prOpts->bPercID)) {
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
*
*/
/* non-residue/gap characters will crash AlnToHMM2(),
therefore sanitise unknown characters, FS, r259 -> r260 */
SanitiseUnknown(prMSeq);
if (iIterationCounter < prOpts->iMaxHMMIterations) {
Log(&rLog, LOG_INFO, "Computing HMM from alignment");
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,
prOpts->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, prOpts->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
* @param[in] rHhalignPara
* FIXME
*
* @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() */
|
