+++ /dev/null
-/* -*- 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: hhalign_wrapper.c 256 2011-06-23 13:57:28Z fabian $
- */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <stdlib.h>
-#include <string.h>
-#include <assert.h>
-#include <ctype.h>
-
-#include "seq.h"
-#include "tree.h"
-#include "progress.h"
-#include "hhalign/general.h"
-#include "hhalign/hhfunc.h"
-#include "hhalign/hhalign.h"
-
-/* up to this level (from leaf) will background HMM info be applied */
-#define APPLY_BG_HMM_UP_TO_TREE_DEPTH 10
-
-#define TIMING 0
-
-#define TRACE 0
-
-/**
- * @brief get rid of unknown residues
- *
- * @note HHalignWrapper can be entered in 2 different ways: (i) all
- * sequences are un-aligned (ii) there are 2 (aligned) profiles. in
- * the un-aligned case (i) the sequences come straight from Squid,
- * that is, they have been sanitised, all non-alphabetic residues
- * have been rendered as X's. In profile mode (ii) one profile may
- * have been produced internally. In that case residues may have
- * been translated back into their 'native' form, that is, they may
- * contain un-sanitised residues. These will cause trouble during
- * alignment
- * FS, r213->214
- */
-void
-SanitiseUnknown(mseq_t *mseq)
-{
-
- int iS; /* iterator for sequence */
- int iR; /* iterator for residue */
- int iLen; /* length of sequence */
- char *pcRes = NULL;
-
-
- for (iS = 0; iS < mseq->nseqs; iS++){
-
- for (pcRes = mseq->seq[iS]; '\0' != *pcRes; pcRes++){
-
- if (isgap(*pcRes)){
- continue;
- }
-
- if (mseq->seqtype==SEQTYPE_PROTEIN) {
- if (NULL == strchr(AMINO_ALPHABET, toupper(*pcRes))) {
- *pcRes = AMINOACID_ANY;
- }
- } else if (mseq->seqtype==SEQTYPE_DNA) {
- if (NULL == strchr(DNA_ALPHABET, toupper(*pcRes))) {
- *pcRes = NUCLEOTIDE_ANY;
- }
- } else if (mseq->seqtype==SEQTYPE_RNA) {
- if (NULL == strchr(RNA_ALPHABET, toupper(*pcRes))) {
- *pcRes = NUCLEOTIDE_ANY;
- }
- }
-
- } /* !EO String */
-
- } /* 0 <= iS < mseq->nseqs */
-
- return;
-
-} /*** end: SanitiseUnknown() ***/
-
-/**
- * @brief translate unknown residues back to ambiguity codes;
- * hhalign translates ambiguity codes (B,Z) into unknown residue (X).
- * we still have the original (un-aligned) residue information,
- * by iterating along the original and aligned sequences we can
- * reconstruct where codes have been changed and restore them
- * to their original value
- *
- * @param[in,out] mseq
- * sequence/profile data, mseq->seq [in,out] is changed to conform
- * with mseq->orig_seq [in]
- *
- */
-void
-TranslateUnknown2Ambiguity(mseq_t *mseq)
-{
-
- int iS; /* iterator for sequence */
- int iR, iRo; /* iterator for residue (original) */
- int iChange, iCase, iAmbi; /* counts how many replacements */
- static int siOffset = 'a' - 'A';
-
- for (iS = 0; iS < mseq->nseqs; iS++){
-
- iR = iRo = 0;
- iChange = iCase = iAmbi = 0;
-
- while(('\0' != mseq->seq[iS][iR]) &&
- ('\0' != mseq->orig_seq[iS][iRo])) {
-
- /* skip gaps in aligned sequences */
- while(isgap(mseq->seq[iS][iR])) {
- iR++;
- } /* was gap in unaligned seq
- * this should probably not happen */
- while(isgap(mseq->orig_seq[iS][iRo])) {
- iRo++;
- } /* was gap in aligned seq */
-
- /* check if we reached the end of the sequence after
- * skipping the gaps */
- if ( ('\0' == mseq->seq[iS][iR]) ||
- ('\0' == mseq->orig_seq[iS][iRo]) ){
- break;
- }
-
-
- if (mseq->seq[iS][iR] != mseq->orig_seq[iS][iRo]){
- /* FIXME: count replacements, discard case changes */
- iChange++;
- if ( (mseq->seq[iS][iR] == mseq->orig_seq[iS][iRo]+siOffset) ||
- (mseq->seq[iS][iR] == mseq->orig_seq[iS][iRo]-siOffset) ){
- iCase++;
- }
- else {
- iAmbi++;
- }
-#if TRACE
- Log(&rLog, LOG_FORCED_DEBUG, "seq=%d, pos=(%d:%d), (%c:%c)",
- iS, iR, iRo,
- mseq->seq[iS][iR], mseq->orig_seq[iS][iRo]);
-#endif
- mseq->seq[iS][iR] = mseq->orig_seq[iS][iRo];
- }
-
- iR++;
- iRo++;
-
- } /* !EO seq */
-
- Log(&rLog, LOG_DEBUG,
- "in seq %d re-translated %d residue codes (%d true, %d case)",
- iS, iChange, iAmbi, iCase);
-
- /* assert that both sequences (un/aligned) have terminated */
- /* skip gaps in aligned sequences */
- while(isgap(mseq->seq[iS][iR])) {
- iR++;
- } /* was gap in unaligned seq
- * this should probably not happen */
- while(isgap(mseq->orig_seq[iS][iRo])) {
- iRo++;
- } /* was gap in aligned seq */
-
-
- if ( ('\0' != mseq->seq[iS][iR]) ||
- ('\0' != mseq->orig_seq[iS][iRo]) ){
-
- Log(&rLog, LOG_FATAL, "inconsistency in un/aligned sequence %d\n>%s\n>%s\n",
- iS, mseq->seq[iS], mseq->orig_seq[iS]);
- }
-
- } /* 0 <= iS < mseq->nseqs */
-
-} /*** end: TranslateUnknown2Ambiguity() ***/
-
-
-/**
- * @brief re-attach leading and trailing gaps to alignment
- *
- * @param[in,out] prMSeq,
- * alignment structure (at this stage there should be no un-aligned sequences)
- * @param[in] iProfProfSeparator
- * gives sizes of input profiles, -1 if no input-profiles but un-aligned sequences
- *
- * @note leading and tailing profile columns
- * that only contain gaps have no effect on the alignment
- * and are removed during the alignment. if they are
- * encountered a warning message is printed to screen.
- * some users like to preserve these gap columns
- * FS, r213->214
- */
-void
-ReAttachLeadingGaps(mseq_t *prMSeq, int iProfProfSeparator)
-{
-
- int i, j;
- int iSize1 = 0; /* #seqs in 1st profile */
- int iSize2 = 0; /* #seqs in 2nd profile */
- int iPPS = iProfProfSeparator;
- int iLeadO1 = 0; /* leading gaps in 1st seq of 1st profile */
- int iLeadO2 = 0; /* leading gaps in 1st seq of 2nd profile */
- int iLeadA1 = 0; /* leading gaps in 1st seq of final alignment */
- int iLeadA2 = 0; /* leading gaps in PPS seq of final alignment */
- int iTrailO1 = 0; /* trailing gaps in 1st seq of 1st profile */
- int iTrailO2 = 0; /* trailing gaps in 1st seq of 2nd profile */
- int iTrailA1 = 0; /* trailing gaps in 1st seq of final alignment */
- int iTrailA2 = 0; /* trailing gaps in PPS seq of final alignment */
- int iLen = 0; /* length of final alignment */
- int iLen1 = 0; /* length of 1st profile */
- int iLen2 = 0; /* length of 2nd profile */
- int iCutHead = 0; /* make up truncation at head */
- int iCutTail = 0; /* make up truncation at tail */
- char *pcIter = NULL;
-
- if (-1 == iProfProfSeparator){
- return;
- }
- else {
- assert(prMSeq->aligned);
- assert(prMSeq->nseqs > iProfProfSeparator);
- }
- iSize1 = iProfProfSeparator;
- iSize2 = prMSeq->nseqs - iProfProfSeparator;
- iLen = strlen(prMSeq->seq[0]);
- iLen1 = strlen(prMSeq->orig_seq[0]);
- iLen2 = strlen(prMSeq->orig_seq[iPPS]);
-
- /* count leading/trailing gaps in 1st sequence of 1st/2nd profile and final alignmant */
- for (iLeadO1 = 0, pcIter = prMSeq->orig_seq[0]; isgap(*pcIter); pcIter++, iLeadO1++);
- for (iLeadO2 = 0, pcIter = prMSeq->orig_seq[iPPS]; isgap(*pcIter); pcIter++, iLeadO2++);
- for (iLeadA1 = 0, pcIter = prMSeq->seq[0]; isgap(*pcIter); pcIter++, iLeadA1++);
- for (iLeadA2 = 0, pcIter = prMSeq->seq[iPPS]; isgap(*pcIter); pcIter++, iLeadA2++);
- for (iTrailO1 = 0, pcIter = &prMSeq->orig_seq[0][iLen1-1]; isgap(*pcIter); pcIter--, iTrailO1++);
- for (iTrailO2 = 0, pcIter = &prMSeq->orig_seq[iPPS][iLen2-1]; isgap(*pcIter); pcIter--, iTrailO2++);
- for (iTrailA1 = 0, pcIter = &prMSeq->seq[0][iLen-1]; isgap(*pcIter); pcIter--, iTrailA1++);
- for (iTrailA2 = 0, pcIter = &prMSeq->seq[iPPS][iLen-1]; isgap(*pcIter); pcIter--, iTrailA2++);
-
-
- /* turn leading/trailing gaps into truncation */
- iLeadO1 = iLeadO1 > iLeadA1 ? iLeadO1-iLeadA1 : 0;
- iLeadO2 = iLeadO2 > iLeadA2 ? iLeadO2-iLeadA2 : 0;
- iTrailO1 = iTrailO1 > iTrailA1 ? iTrailO1-iTrailA1 : 0;
- iTrailO2 = iTrailO2 > iTrailA2 ? iTrailO2-iTrailA2 : 0;
- iCutHead = iLeadO1 > iLeadO2 ? iLeadO1 : iLeadO2;
- iCutTail = iTrailO1 > iTrailO2 ? iTrailO1 : iTrailO2;
-
- /* re-allocate and shift memory */
- if ( (iCutHead > 0) || (iCutTail > 0) ){ /* skip if no re-attachment, FS, r244 -> r245 */
- for (i = 0; i < prMSeq->nseqs; i++){
-
- CKREALLOC(prMSeq->seq[i], iLen+iCutHead+iCutTail+2);
- if (iCutHead > 0){ /* skip if no re-attachment, FS, r244 -> r245 */
- memmove(prMSeq->seq[i]+iCutHead, prMSeq->seq[i], iLen);
- }
- for (j = 0; j < iCutHead; j++){
- prMSeq->seq[i][j] = '-';
- }
- for (j = iLen+iCutHead; j < iLen+iCutHead+iCutTail; j++){
- prMSeq->seq[i][j] = '-';
- }
- prMSeq->seq[i][j] = '\0';
- }
- } /* (iCutHead > 0) || (iCutTail > 0) */
-
-} /*** end: ReAttachLeadingGaps() ***/
-
-/**
- * @brief reallocate enough memory for alignment and
- * attach sequence pointers to profiles
- *
- * @param[in,out] mseq
- * sequence/profile data, increase memory for sequences in profiles
- * @param[out] ppcProfile1
- * pointers to sequencese in 1st profile
- * @param[out] ppcProfile2
- * pointers to sequencese in 2nd profile
- * @param[out] pdWeightsL
- * weights (normalised to 1.0) for sequences in left profile
- * @param[out] pdWeightsR
- * weights (normalised to 1.0) for sequences in right profile
- * @param[in] pdSeqWeights
- * weights for _all_ sequences in alignment
- * @param[in] iLeafCountL
- * number of sequences in 1st profile
- * @param[in] piLeafListL
- * array of integer IDs of sequences in 1st profile
- * @param[in] iLeafCountR
- * number of sequences in 2nd profile
- * @param[in] piLeafListR
- * array of integer IDs of sequences in 2nd profile
- *
- */
-void
-PrepareAlignment(mseq_t *mseq, char **ppcProfile1, char **ppcProfile2,
- double *pdWeightsL, double *pdWeightsR, double *pdSeqWeights,
- int iLeafCountL, int *piLeafListL,
- int iLeafCountR, int *piLeafListR)
-{
-
- int iLenL = 0; /* length of 1st profile */
- int iLenR = 0; /* length of 2nd profile */
- int iMaxLen = 0; /* maximum possible length of alignment */
- int i; /* aux */
- double dWeight = 0.00;
- double dWeightInv = 0.00;
-
- assert(NULL!=mseq);
- assert(NULL!=ppcProfile1);
- assert(NULL!=ppcProfile2);
- assert(NULL!=piLeafListL);
- assert(NULL!=piLeafListR);
-
- /* get length of profiles,
- * in a profile all sequences should have same length so only look at 1st
- */
- iLenL = strlen(mseq->seq[piLeafListL[0]]);
- iLenR = strlen(mseq->seq[piLeafListR[0]]);
- iMaxLen = iLenL + iLenR + 1;
-
- /* reallocate enough memory for sequences in alignment (for adding
- * gaps)
- */
- for (i = 0; i < iLeafCountL; i++){
- mseq->seq[piLeafListL[i]] =
- CKREALLOC(mseq->seq[piLeafListL[i]], iMaxLen);
- }
- for (i = 0; i < iLeafCountR; i++){
- mseq->seq[piLeafListR[i]] =
- CKREALLOC(mseq->seq[piLeafListR[i]], iMaxLen);
- }
-
- /* attach sequences to profiles
- */
- for (i = 0; i < iLeafCountL; i++){
- ppcProfile1[i] = mseq->seq[piLeafListL[i]];
- }
- ppcProfile1[i] = NULL;
-
- for (i = 0; i < iLeafCountR; i++){
- ppcProfile2[i] = mseq->seq[piLeafListR[i]];
- }
- ppcProfile2[i] = NULL;
-
-
- /* remove terminal 'X' for single sequences:
- * it is a quirk of hhalign() to delete 2 individual sequences
- * if 2 terminal X's meet during alignment,
- * just replace (one of) them.
- * this can be undone at the end.
- * profiles -consisting of more than 1 sequence-
- * appear to be all-right.
- * there seems to be no problem with B's and Z's
- */
- if ( (1 == iLeafCountL) && (1 == iLeafCountR) ){
-
- if ( ('X' == ppcProfile1[0][0]) && ('X' == ppcProfile2[0][0]) ){
-#define NOTX 'N'
- ppcProfile1[0][0] = NOTX; /* FIXME: arbitrary assignment */
- ppcProfile2[0][0] = NOTX; /* FIXME: arbitrary assignment */
- }
- if ( ('X' == ppcProfile1[0][iLenL-1]) && ('X' == ppcProfile2[0][iLenR-1]) ){
- ppcProfile1[0][iLenL-1] = NOTX; /* FIXME: arbitrary assignment */
- ppcProfile2[0][iLenR-1] = NOTX; /* FIXME: arbitrary assignment */
- }
- }
-
- /* obtain sequence weights
- */
- if (NULL != pdSeqWeights){
-
- dWeight = 0.00;
- for (i = 0; i < iLeafCountL; i++){
- register double dAux = pdSeqWeights[piLeafListL[i]];
-#ifndef NDEBUG
- if (dAux <= 0.00){
- Log(&rLog, LOG_DEBUG, "seq-weight %d = %f", piLeafListL[i], dAux);
- }
-#endif
- pdWeightsL[i] = dAux;
- dWeight += dAux;
- } /* 0 <= i < iLeafCountL */
- dWeightInv = 1.00 / dWeight;
- for (i = 0; i < iLeafCountL; i++){
- pdWeightsL[i] *= dWeightInv;
- }
-
- dWeight = 0.00;
- for (i = 0; i < iLeafCountR; i++){
- register double dAux = pdSeqWeights[piLeafListR[i]];
-#ifndef NDEBUG
- if (dAux <= 0.00){
- Log(&rLog, LOG_DEBUG, "seq-weight %d = %f", piLeafListR[i], dAux);
- }
-#endif
- pdWeightsR[i] = dAux;
- dWeight += dAux;
- } /* 0 <= i < iLeafCountL */
- dWeightInv = 1.00 / dWeight;
- for (i = 0; i < iLeafCountR; i++){
- pdWeightsR[i] *= dWeightInv;
- }
- } /* (NULL != pdSeqWeights) */
- else {
- pdWeightsL[0] = pdWeightsR[0] = -1.00;
- }
-
-#if TRACE
- for (i = 0; i < iLeafCountL; i++){
- Log(&rLog, LOG_FORCED_DEBUG, "ppcProfile1[%d/%d] pointing to mseq %d = %s",
- i, iLeafCountR, piLeafListL[i], ppcProfile1[i]);
- }
- for (i = 0; i < iLeafCountR; i++){
- Log(&rLog, LOG_FORCED_DEBUG, "ppcProfile2[%d/%d] pointing to mseq %d = %s",
- i, iLeafCountR, piLeafListR[i], ppcProfile2[i]);
- }
-#endif
-
- return;
-
-} /*** end: PrepareAlignment() ***/
-
-
-/**
- * @brief wrapper for hhalign. This is a frontend function to
- * the ported hhalign code.
- *
- * @param[in,out] prMSeq
- * holds the unaligned sequences [in] and the final alignment [out]
- * @param[in] piOrderLR
- * holds order in which sequences/profiles are to be aligned,
- * even elements specify left nodes, odd elements right nodes,
- * if even and odd are same then it is a leaf
- * @param[in] pdSeqWeights
- * Weight per sequence. No weights used if NULL
- * @param[in] iNodeCount
- * number of nodes in tree, piOrderLR has 2*iNodeCount elements
- * @param[in] prHMMList
- * List of background HMMs (transition/emission probabilities)
- * @param[in] iHMMCount
- * Number of input background HMMs
- * @param[in] iProfProfSeparator
- * Gives the number of sequences in the first profile, if in
- * profile/profile alignment mode (iNodeCount==3). That assumes mseqs
- * holds the sequences of profile 1 and profile 2.
- * @param[in] rHhalignPara
- * various parameters read from commandline
- *
- * @return score of the alignment FIXME what is this?
- *
- * @note complex function. could use some simplification, more and
- * documentation and a struct'uring of piOrderLR
- *
- * @note HHalignWrapper can be entered in 2 different ways:
- * (i) all sequences are un-aligned (ii) there are 2 (aligned) profiles.
- * in the un-aligned case (i) the sequences come straight from Squid,
- * that is, they have been sanitised, all non-alphabetic residues
- * have been rendered as X's. In profile mode (ii) one profile may
- * have been produced internally. In that case residues may have
- * been translated back into their 'native' form, that is, they
- * may contain un-sanitised residues. These will cause trouble
- * during alignment
- *
- * @note: introduced argument hhalign_para rHhalignPara; FS, r240 -> r241
- * @note: if hhalign() fails then try with Viterbi by setting MAC-RAM=0; FS, r241 -> r243
- */
-
-double
-HHalignWrapper(mseq_t *prMSeq, int *piOrderLR,
- double *pdSeqWeights, int iNodeCount,
- hmm_light *prHMMList, int iHMMCount,
- int iProfProfSeparator, hhalign_para rHhalignPara)
-{
- int iN; /* node iterator */
- int *piLeafCount = NULL; /* number of leaves beneath a certain node */
- int **ppiLeafList = NULL; /* list of leaves beneath a certain node */
- char **ppcProfile1 = NULL; /* pointer to sequences in profile */
- char **ppcProfile2 = NULL; /* pointer to sequences in profile */
- char *pcReprsnt1 = NULL; /* representative of HMM aligned to left */
- char *pcReprsnt2 = NULL; /* representative of HMM aligned to right */
- char **ppcReprsnt1 = &pcReprsnt1; /* representative of HMM aligned to L */
- char **ppcReprsnt2 = &pcReprsnt2; /* representative of HMM aligned to R */
- char *pcConsens1 = NULL; /* copy of left sequence */
- char *pcConsens2 = NULL; /* copy of right sequence */
- char **ppcCopy1 = /*&pcCopy1*/NULL; /* copy of left sequences */
- char **ppcCopy2 = /*&pcCopy2*/NULL; /* copy of right sequences */
- double *pdScores = NULL; /* alignment scores (seq/HMM) */
- double dScore = 0.0; /* alignment score (seq/HMM) */
- int iAux_FS = 0;
- char zcAux[10000] = {0};
- char zcError[10000] = {0};
- int i; /* aux */
- progress_t *prProgress;
- int iAlnLen; /* alignment length */
- double *pdWeightsL = NULL; /* sequence weights of left profile */
- double *pdWeightsR = NULL; /* sequence weights of right profile */
- int iMergeNodeCounter = 0;
- hmm_light *prHMM = NULL;
- bool bPrintCR = (rLog.iLogLevelEnabled<=LOG_VERBOSE) ? FALSE : TRUE;
-#if TIMING
- char zcStopwatchMsg[1024];
- Stopwatch_t *stopwatch = StopwatchCreate();
- StopwatchZero(stopwatch);
- StopwatchStart(stopwatch);
-#endif
-
-
- if (NULL != prHMMList) {
- if (iHMMCount>1) {
- Log(&rLog, LOG_WARN, "FIXME: Using only first of %u HMMs (needs implementation)", iHMMCount);
- }
- prHMM = &(prHMMList[0]);
- } else {
- /* FIXME: prHMM not allowed to be NULL and therefore pseudo allocated here */
- prHMM = (hmm_light *) CKCALLOC(1, sizeof(hmm_light));
- }
-
- assert(NULL!=prMSeq);
-
- if (NULL==piOrderLR) {
- assert(3==iNodeCount);
- }
- SanitiseUnknown(prMSeq);
-
- /* hhalign was not made for DNA/RNA. So warn if sequences are not
- * protein
- */
- if (SEQTYPE_PROTEIN != prMSeq->seqtype) {
- Log(&rLog, LOG_WARN, "Sequence type is %s. HHalign has only been tested on protein.",
- SeqTypeToStr(prMSeq->seqtype));
- }
-
- /* hhalign produces funny results if sequences contain gaps, so
- * dealign. Only way to use alignment info is to use it as a
- * background HMM
- */
- if (TRUE == prMSeq->aligned) {
- Log(&rLog, LOG_DEBUG, "Dealigning aligned sequences (inside %s)", __FUNCTION__);
- DealignMSeq(prMSeq);
- }
-
-#if TRACE
- Log(&rLog, LOG_FORCED_DEBUG, "iNodeCount = %d", iNodeCount);
-#endif
-
-
- /* allocate top-level memory for leaf tracking arrays and profiles,
- * and sequence weights*/
- piLeafCount = CKCALLOC(iNodeCount, sizeof(int));
- ppiLeafList = CKCALLOC(iNodeCount, sizeof(int *));
- ppcProfile1 = CKCALLOC(prMSeq->nseqs*2-1, sizeof(char *));
- ppcProfile2 = CKCALLOC(prMSeq->nseqs*2-1, sizeof(char *));
- pdScores = CKCALLOC(iNodeCount, sizeof(double));
- pdWeightsL = CKCALLOC(iNodeCount, sizeof(double));
- pdWeightsR = CKCALLOC(iNodeCount, sizeof(double));
-
- NewProgress(&prProgress, LogGetFP(&rLog, LOG_INFO),
- "Progressive alignment progress", bPrintCR);
-
-
- /* Profile-profile alignment? Then setup piLeafCount and
- * piLeafList here. FIXME this is just an awful haaaack
- */
- if (iNodeCount==3 && NULL==piOrderLR) {
- int iSizeProf1 = iProfProfSeparator;
- int iSizeProf2 = prMSeq->nseqs - iProfProfSeparator;
-
- piLeafCount[0] = iSizeProf1;
- ppiLeafList[0] = (int *)CKMALLOC(iSizeProf1 * sizeof(int));
- for (i=0;i<iSizeProf1;i++) {
- ppiLeafList[0][i] = i;
- }
- piLeafCount[1] = iSizeProf2;
- ppiLeafList[1] = (int *)CKMALLOC(iSizeProf2 * sizeof(int));
- for (i=0;i<iSizeProf2;i++) {
- ppiLeafList[1][i] = i+iSizeProf1;
- }
-
- /* awful hack inside an awful hack: we had to setup piLeafCount
- * and piLeafList outside the node iteration. this which is
- * normally done at leaf level inside the node iteration. to
- * avoid overwriting the already setup vars set...
- */
- iNodeCount=1;
-
- piOrderLR = (int *)CKCALLOC(DIFF_NODE, sizeof(int));
- piOrderLR[LEFT_NODE] = 0;
- piOrderLR[RGHT_NODE] = 1;
- piOrderLR[PRNT_NODE] = 2;
- }
-
-
-
- iMergeNodeCounter = 0;
- for (iN = 0; iN < iNodeCount; iN++){
-
- register int riAux = DIFF_NODE * iN;
-
- /*LOG_DEBUG("node %d ", iN);*/
-
- if (piOrderLR[riAux+LEFT_NODE] == piOrderLR[riAux+RGHT_NODE]){
-
- register int riLeaf = piOrderLR[riAux+LEFT_NODE];
-#if TRACE
- if (NULL == pdSeqWeights) {
- Log(&rLog, LOG_FORCED_DEBUG, "node %d is a leaf with entry %d (seq %s)",
- iN, riLeaf, prMSeq->sqinfo[riLeaf].name);
- } else {
- Log(&rLog, LOG_FORCED_DEBUG, "node %d is a leaf with entry %d (seq %s) and weight %f",
- iN, riLeaf, prMSeq->sqinfo[riLeaf].name, pdSeqWeights[riLeaf]);
- }
-#endif
- /* left/right entry same, this is a leaf */
- piLeafCount[piOrderLR[riAux+PRNT_NODE]] = 1; /* number of leaves is '1' */
- ppiLeafList[piOrderLR[riAux+PRNT_NODE]] = (int *)CKMALLOC(1 * sizeof(int));
- ppiLeafList[piOrderLR[riAux+PRNT_NODE]][0] = riLeaf;
-
- } /* was a leaf */
- else {
- int iL, iR, iP; /* ID of left/right nodes, parent */
- int i, j; /* aux */
-
- Log(&rLog, LOG_DEBUG,
- "merge profiles at node %d", iN, piOrderLR[riAux]);
-
- /* iNodeCount - prMSeq->nseqs = total # of merge-nodes
- * unless in profile/profile alignment mode
- */
- if (1 == iNodeCount) {
- ProgressLog(prProgress, ++iMergeNodeCounter,
- 1, FALSE);
- } else {
- ProgressLog(prProgress, ++iMergeNodeCounter,
- iNodeCount - prMSeq->nseqs, FALSE);
- }
-
- /* left/right entry are not same, this is a merge node */
- iL = piOrderLR[riAux+LEFT_NODE];
- iR = piOrderLR[riAux+RGHT_NODE];
- iP = piOrderLR[riAux+PRNT_NODE];
- piLeafCount[iP] = piLeafCount[iL] + piLeafCount[iR];
- ppiLeafList[iP] = (int *)CKMALLOC(piLeafCount[iP] * sizeof(int));
-
- for (i = j = 0; i < piLeafCount[iL]; i++, j++){
- ppiLeafList[iP][j] = ppiLeafList[iL][i];
- }
- for (i = 0; i < piLeafCount[iR]; i++, j++){
- ppiLeafList[iP][j] = ppiLeafList[iR][i];
- }
-
- /* prepare simulation arena:
- * - make sure enough memory in sequences
- * - attach sequence pointers to profiles
- */
- /* idea: switch template and query according to nseq? */
-
- PrepareAlignment(prMSeq, ppcProfile1, ppcProfile2,
- pdWeightsL, pdWeightsR, pdSeqWeights,
- piLeafCount[iL], ppiLeafList[iL],
- piLeafCount[iR], ppiLeafList[iR]);
- if (rLog.iLogLevelEnabled <= LOG_DEBUG){
- int i;
- FILE *fp = LogGetFP(&rLog, LOG_DEBUG);
- Log(&rLog, LOG_DEBUG, "merging profiles %d & %d", iL, iR);
- for (i = 0; i < piLeafCount[iL]; i++){
- fprintf(fp, "L/#=%3d (ID=%3d, w=%f): %s\n",
- i, ppiLeafList[iL][i], pdWeightsL[i], ppcProfile1[i]);
- }
- for (i = 0; i < piLeafCount[iR]; i++){
- fprintf(fp, "R/#=%3d (ID=%3d, w=%f): %s\n",
- i, ppiLeafList[iR][i], pdWeightsR[i], ppcProfile2[i]);
- }
- }
-
-
-
- /* align individual sequences to HMM;
- * - use representative sequence to get gapping
- * - create copies of both, individual/representative sequences
- * as we don't want to introduce gaps into original
- *
- * FIXME: representative sequence is crutch, should use
- * full HMM but that does not seem to work at all
- * -- try harder! Fail better!
- */
- if ( (piLeafCount[iL] <= APPLY_BG_HMM_UP_TO_TREE_DEPTH) && (0 != prHMM->L) ){
- int i, j;
- pcReprsnt1 = CKCALLOC(prHMM->L+strlen(ppcProfile1[0])+1, sizeof(char));
- for (i = 0; i < prHMM->L; i++){
- pcReprsnt1[i] = prHMM->seq[prHMM->ncons][i+1];
- }
- ppcCopy1 = CKCALLOC(piLeafCount[iL], sizeof(char *));
- for (j = 0; j < piLeafCount[iL]; j++){
- ppcCopy1[j] = CKCALLOC(prHMM->L+strlen(ppcProfile1[0])+1, sizeof(char));
- for (i = 0; i < (int) strlen(ppcProfile1[0]); i++){
- ppcCopy1[j][i] = ppcProfile1[j][i];
- }
- }
-
- {
- /* the size of the elements in the forward/backward matrices
- depends very much on the lengths of the profiles _and_
- in which position (1st/2nd) the longer/shorter profile/HMM is.
- the matrix elements can easily exceed the size of a (long?) double
- if the longer profile/HMM is associated with the query (q) and the
- shorter with the target (t).
- FIXME: however, pseudo-count adding may also depend on position,
- this is only really tested for the HMM being in the 1st position (q)
- MUST TEST THIS MORE THOROUGHLY
-
- this switch appears to be most easily (although unelegantly)
- effected here. Don't want to do it (upstairs) in PrepareAlignment()
- because it might jumble up the nodes. Don't want to do it in hhalign()
- either because ppcProfile1/2 and q/t may be used independently.
- FS, r236 -> r237
- */
- int iLenA = strlen(ppcCopy1[0]);
- int iLenH = prHMM->L;
- int iHHret = 0;
-
- if (iLenH < iLenA){
- iHHret = hhalign(ppcReprsnt1, 0/* only one representative seq*/, NULL,
- ppcCopy1, piLeafCount[iL], pdWeightsL,
- &dScore, prHMM,
- NULL, NULL, NULL, NULL,
- rHhalignPara,
- iAux_FS++, /* DEBUG ARGUMENT */
- zcAux, zcError);
- }
- else {
- iHHret = hhalign(ppcCopy1, piLeafCount[iL], pdWeightsL,
- ppcReprsnt1, 0/* only one representative seq*/, NULL,
- &dScore, prHMM,
- NULL, NULL, NULL, NULL,
- rHhalignPara,
- iAux_FS++, /* DEBUG ARGUMENT */
- zcAux, zcError);
- }
- if (0 != iHHret){ /* FS, r255 -> */
- fprintf(stderr, "%s:%s:%d: (not essential) HMM pre-alignment failed, error %d, \n"
- "\t#=%d (len=%d), lead-seq=%s, len(HMM)=%d\tCARRY ON REGARDLESS\n",
- __FUNCTION__, __FILE__, __LINE__, iHHret,
- piLeafCount[iL], (int)strlen(ppcCopy1[0]), prMSeq->sqinfo[ppiLeafList[iL][0]].name, (int)strlen(ppcReprsnt1[0]));
- }
- }
- pdScores[ppiLeafList[iL][0]] = dScore;
-#if 0
- printf("score: %f\nL: %s\nH: %s\n",
- dScore, ppcCopy1[0], ppcReprsnt1[0]);
-#endif
- /* assemble 'consensus';
- * this is not a real consensus, it is more a gap indicator,
- * for each position it consists of residues/gaps in the 1st sequences,
- * or a residue (if any) of the other sequences.
- * it only contains a gap if all sequences of the profile
- * have a gap at this position
- */
- pcConsens1 = CKCALLOC(prHMM->L+strlen(ppcProfile1[0])+1, sizeof(char));
- for (i = 0; i < prHMM->L; i++){
- for (j = 0, pcConsens1[i] = '-'; (j < piLeafCount[iL]) && ('-' == pcConsens1[i]); j++){
- pcConsens1[i] = ppcCopy1[j][i];
- }
- }
-#if 0
- for (j = 0; (j < piLeafCount[iL]); j++){
- printf("L%d:%s\n", j, ppcCopy1[j]);
- }
- printf("LC:%s\n", pcConsens1);
-#endif
- } /* ( (1 == piLeafCount[iL]) && (0 != prHMM->L) ) */
-
- if ( (piLeafCount[iR] <= APPLY_BG_HMM_UP_TO_TREE_DEPTH) && (0 != prHMM->L) ){
- int i, j;
-
- pcReprsnt2 = CKCALLOC(prHMM->L+strlen(ppcProfile2[0])+1, sizeof(char));
- for (i = 0; i < prHMM->L; i++){
- pcReprsnt2[i] = prHMM->seq[prHMM->ncons][i+1];
- }
- ppcCopy2 = CKCALLOC(piLeafCount[iR], sizeof(char *));
- for (j = 0; j < piLeafCount[iR]; j++){
- ppcCopy2[j] = CKCALLOC(prHMM->L+strlen(ppcProfile2[0])+1, sizeof(char));
- for (i = 0; i < (int) strlen(ppcProfile2[0]); i++){
- ppcCopy2[j][i] = ppcProfile2[j][i];
- }
- }
-
- {
- /* the size of the elements in the forward/backward matrices
- depends very much on the lengths of the profiles _and_
- in which position (1st/2nd) the longer/shorter profile/HMM is.
- the matrix elements can easily exceed the size of a (long?) double
- if the longer profile/HMM is associated with the query (q) and the
- shorter with the target (t).
- FIXME: however, pseudo-count adding may also depend on position,
- this is only really tested for the HMM being in the 1st position (q)
- MUST TEST THIS MORE THOROUGHLY
-
- this switch appears to be most easily (although unelegantly)
- effected here. Don't want to do it (upstairs) in PrepareAlignment()
- because it might jumble up the nodes. Don't want to do it in hhalign()
- either because ppcProfile1/2 and q/t may be used independently.
- FS, r236 -> r237
- */
- int iLenA = strlen(ppcCopy2[0]);
- int iLenH = prHMM->L;
- int iHHret = 0;
-
- if (iLenH < iLenA){
- iHHret = hhalign(ppcReprsnt2, 0/* only one representative seq */, NULL,
- ppcCopy2, piLeafCount[iR], pdWeightsR,
- &dScore, prHMM,
- NULL, NULL, NULL, NULL,
- rHhalignPara,
- iAux_FS++, /* DEBUG ARGUMENT */
- zcAux, zcError);
- }
- else {
- iHHret = hhalign(ppcCopy2, piLeafCount[iR], pdWeightsR,
- ppcReprsnt2, 0/* only one representative seq */, NULL,
- &dScore, prHMM,
- NULL, NULL, NULL, NULL,
- rHhalignPara,
- iAux_FS++, /* DEBUG ARGUMENT */
- zcAux, zcError);
- }
- if (0 != iHHret){ /* FS, r255 -> */
- fprintf(stderr, "%s:%s:%d: (not essential) HMM pre-alignment failed, error %d, \n"
- "\t#=%d (len=%d), lead-seq=%s, len(HMM)=%d\tCARRY ON REGARDLESS\n",
- __FUNCTION__, __FILE__, __LINE__, iHHret,
- piLeafCount[iR], (int)strlen(ppcCopy2[0]), prMSeq->sqinfo[ppiLeafList[iR][0]].name, (int)strlen(ppcReprsnt2[0]));
- }
- }
- pdScores[ppiLeafList[iR][0]] = dScore;
-#if 0
- printf("H: %s\nR: %s\nscore: %f\n",
- ppcReprsnt2[0], ppcCopy2[0], dScore);
-#endif
- /* assemble 'consensus';
- * this is not a real consensus, it is more a gap indicator,
- * for each position it consists of residues/gaps in the 1st sequences,
- * or a residue (if any) of the other sequences.
- * it only contains a gap if all sequences of the profile
- * have a gap at this position
- */
- pcConsens2 = CKCALLOC(prHMM->L+strlen(ppcProfile2[0])+1, sizeof(char));
- for (i = 0; i < prHMM->L; i++){
- for (j = 0, pcConsens2[i] = '-'; (j < piLeafCount[iR]) && ('-' == pcConsens2[i]); j++){
- pcConsens2[i] = ppcCopy2[j][i];
- }
- }
-#if 0
- for (j = 0; (j < piLeafCount[iR]); j++){
- printf("R%d:%s\n", j, ppcCopy2[j]);
- }
- printf("RC:%s\n", pcConsens2);
-#endif
- } /* ( (1 == piLeafCount[iR]) && (0 != prHMM->L) ) */
-
-
-
- /* do alignment here (before free)
- */
- {
- /* the size of the elements in the forward/backward matrices
- depends very much on the lengths of the profiles _and_
- in which position (1st/2nd) the longer/shorter profile is.
- the matrix elements can easily exceed the size of a (long?) double
- if the longer profile is associated with the query (q) and the
- shorter with the target (t).
- this switch appears to be most easily (although unelegantly)
- effected here. Don't want to do it (upstairs) in PrepareAlignment()
- because it might jumble up the nodes. Don't want to do it in hhalign()
- either because ppcProfile1/2 and q/t may be used independently.
- FS, r228 -> 229
- */
- int iLen1 = strlen(ppcProfile1[0]);
- int iLen2 = strlen(ppcProfile2[0]);
- /* potential problem with empty profiles, FS, r249 -> r250 */
- if ( (0 == iLen1) || (0 == iLen2) ){
- Log(&rLog, LOG_FATAL, "strlen(prof1)=%d, strlen(prof2)=%d -- nothing to align\n",
- iLen1, iLen2);
- }
-
-
- if (iLen1 < iLen2){
- int iHHret = 0;
- int iOldMacRam = rHhalignPara.iMacRamMB;
- iHHret = hhalign(ppcProfile1, piLeafCount[iL], pdWeightsL,
- ppcProfile2, piLeafCount[iR], pdWeightsR,
- &dScore, prHMM,
- pcConsens1, pcReprsnt1,
- pcConsens2, pcReprsnt2,
- rHhalignPara,
- iAux_FS++, /* DEBUG ARGUMENT */
- zcAux, zcError);
-
- if (RETURN_OK != iHHret){ /* FS, r241 -> */
- fprintf(stderr,
- "%s:%s:%d: problem in alignment (profile sizes: %d + %d) (%s + %s), forcing Viterbi\n"
- "\thh-error-code=%d (mac-ram=%d)\n",
- __FUNCTION__, __FILE__, __LINE__, piLeafCount[iL], piLeafCount[iR],
- prMSeq->sqinfo[ppiLeafList[iL][0]].name, prMSeq->sqinfo[ppiLeafList[iR][0]].name,
- iHHret, rHhalignPara.iMacRamMB);
- /* at this stage hhalign() has failed,
- the only thing we can do (easily) is to re-run it in Viterbi mode,
- for this set MAC-RAM=0, set it back to its original value after 2nd try.
- FS, r241 -> r243 */
- if (RETURN_FROM_MAC == iHHret){
- /* Note: the default way to run hhalign() is to initially select MAC
- by giving it all the memory it needs. MAC may fail due to overflow (repeats?).
- alternatively, the problem may be (genuinely) too big for MAC.
- in thses cases it is legitimate to switch to Viterbi.
- However, selecting Viterbi from the outset is an abuse (abomination!),
- should this 1st invocation of Viterbi fail, then we (FS) will overrule
- the user and hammer the system with a massive memory request.
- (Jos 2:19) If anyone goes outside your house into the street,
- his blood will be on his own head; we will not be responsible. FS, r246 -> r247 */
- rHhalignPara.iMacRamMB = 0;
- }
- else {
- rHhalignPara.iMacRamMB = REALLY_BIG_MEMORY_MB;
- }
-
- iHHret = hhalign(ppcProfile1, piLeafCount[iL], pdWeightsL,
- ppcProfile2, piLeafCount[iR], pdWeightsR,
- &dScore, prHMM,
- pcConsens1, pcReprsnt1,
- pcConsens2, pcReprsnt2,
- rHhalignPara,
- iAux_FS++, /* DEBUG ARGUMENT */
- zcAux, zcError);
- if (RETURN_OK != iHHret){ /* at this stage hhalign() has failed twice,
- 1st time MAC, 2nd time Viterbi, don't know what to do else. FS, r241 -> r243 */
- fprintf(stderr, "%s:%s:%d: problem in alignment, even Viterbi did not work\n"
- "\thh-error-code=%d (mac-ram=%d)\n",
- __FUNCTION__, __FILE__, __LINE__, iHHret, rHhalignPara.iMacRamMB);
- Log(&rLog, LOG_FATAL, "could not perform alignment -- bailing out\n");
- }
- else {
- fprintf(stderr, "%s:%s:%d: 2nd attempt worked", __FUNCTION__, __FILE__, __LINE__);
- }
- rHhalignPara.iMacRamMB = iOldMacRam;
- } /* 1st invocation failed */
-
- } /* 1st profile was shorter than 2nd */
- else {
- int iHHret = 0;
- int iOldMacRam = rHhalignPara.iMacRamMB;
- iHHret = hhalign(ppcProfile2, piLeafCount[iR], pdWeightsR,
- ppcProfile1, piLeafCount[iL], pdWeightsL,
- &dScore, prHMM,
- pcConsens2, pcReprsnt2,
- pcConsens1, pcReprsnt1,
- rHhalignPara,
- iAux_FS++, /* DEBUG ARGUMENT */
- zcAux, zcError);
-
- if (RETURN_OK != iHHret){ /* FS, r241 -> r243 */
- fprintf(stderr,
- "%s:%s:%d: problem in alignment (profile sizes: %d + %d) (%s + %s), forcing Viterbi\n"
- "\thh-error-code=%d (mac-ram=%d)\n",
- __FUNCTION__, __FILE__, __LINE__, piLeafCount[iL], piLeafCount[iR],
- prMSeq->sqinfo[ppiLeafList[iL][0]].name, prMSeq->sqinfo[ppiLeafList[iR][0]].name,
- iHHret, rHhalignPara.iMacRamMB);
- /* at this stage hhalign() has failed,
- the only thing we can do (easily) is to re-run it in Viterbi mode,
- for this set MAC-RAM=0, set it back to its original value after 2nd try.
- FS, r241 -> r243 */
- if (RETURN_FROM_MAC == iHHret){
- /* see above */
- rHhalignPara.iMacRamMB = 0;
- }
- else {
- rHhalignPara.iMacRamMB = REALLY_BIG_MEMORY_MB;
- }
-
- iHHret = hhalign(ppcProfile2, piLeafCount[iR], pdWeightsR,
- ppcProfile1, piLeafCount[iL], pdWeightsL,
- &dScore, prHMM,
- pcConsens2, pcReprsnt2,
- pcConsens1, pcReprsnt1,
- rHhalignPara,
- iAux_FS++, /* DEBUG ARGUMENT */
- zcAux, zcError);
- if (RETURN_OK != iHHret){ /* at this stage hhalign() has failed twice,
- 1st time MAC, 2nd time Viterbi, don't know what to do else. FS, r241 -> r243 */
- fprintf(stderr, "%s:%s:%d: problem in alignment, even Viterbi did not work\n"
- "\thh-error-code=%d (mac-ram=%d)\n",
- __FUNCTION__, __FILE__, __LINE__, iHHret, rHhalignPara.iMacRamMB);
- Log(&rLog, LOG_FATAL, "could not perform alignment -- bailing out\n");
- }
- else {
- fprintf(stderr, "%s:%s:%d: 2nd attempt worked", __FUNCTION__, __FILE__, __LINE__);
- }
- rHhalignPara.iMacRamMB = iOldMacRam;
- } /* 1st invocation failed */
-
- } /* 2nd profile was shorter than 1st */
-
- }
- /* free left/right node lists,
- * after alignment left/right profiles no longer needed
- */
- if (NULL != ppcCopy1){
- int i;
- for (i = 0; i < piLeafCount[iL]; i++){
- CKFREE(ppcCopy1[i]);
- }
- CKFREE(ppcCopy1);
- CKFREE(pcReprsnt1);
- CKFREE(pcConsens1);
- }
- if (NULL != ppcCopy2){
- int i;
- for (i = 0; i < piLeafCount[iR]; i++){
- CKFREE(ppcCopy2[i]);
- }
- CKFREE(ppcCopy2);
- CKFREE(pcReprsnt2);
- CKFREE(pcConsens2);
- }
- ppiLeafList[iL] = CKFREE(ppiLeafList[iL]);
- ppiLeafList[iR] = CKFREE(ppiLeafList[iR]);
- piLeafCount[iL] = piLeafCount[iR] = 0;
-
- } /* was a merge node */
-
- if (rLog.iLogLevelEnabled <= LOG_DEBUG){
- int i, j;
- FILE *fp = LogGetFP(&rLog, LOG_DEBUG);
- for (i = 0; i < iNodeCount; i++){
- if (0 == piLeafCount[i]){
- continue;
- }
- fprintf(fp, "node %3d, #leaves=%d:\t", i, piLeafCount[i]);
- for (j = 0; ppiLeafList && (j < piLeafCount[i]); j++){
- fprintf(fp, "%d,", ppiLeafList[i][j]);
- }
- fprintf(fp, "\n");
- }
- }
-
-
- } /* 0 <= iN < iNodeCount */
- ProgressDone(prProgress);
-
-
- /* check length and set length info
- */
- iAlnLen = strlen(prMSeq->seq[0]);
- for (i=0; i<prMSeq->nseqs; i++) {
-#if 0
- Log(&rLog, LOG_FORCED_DEBUG, "seq no %d: name %s; len %d; %s",
- i, prMSeq->sqinfo[i].name, strlen(prMSeq->seq[i]), prMSeq->seq[i]);
-#endif
-
-#ifndef NDEBUG
- assert(iAlnLen == strlen(prMSeq->seq[i]));
-#endif
- prMSeq->sqinfo[i].len = iAlnLen;
- }
- prMSeq->aligned = TRUE;
-
-
- if (rLog.iLogLevelEnabled <= LOG_DEBUG){
- if (0 != prHMM->L){
- int i;
- Log(&rLog, LOG_DEBUG, "Alignment scores with HMM:");
- for (i = 0; /*pdScores[i] > 0.0*/i < prMSeq->nseqs; i++){
- Log(&rLog, LOG_DEBUG, "%2d:\t%f\n", i, pdScores[i]);
- }
- }
- }
-
-
- /** translate back ambiguity residues
- * hhalign translates ambiguity codes (B,Z) into unknown residues (X).
- * as we still have the original input we can substitute them back
- */
- TranslateUnknown2Ambiguity(prMSeq);
- ReAttachLeadingGaps(prMSeq, iProfProfSeparator);
-
- if (NULL == prHMMList){
- CKFREE(prHMM);
- }
- CKFREE(ppcProfile2);
- CKFREE(ppcProfile1);
- CKFREE(ppiLeafList[piOrderLR[DIFF_NODE*(iNodeCount-1)+PRNT_NODE]]);
- CKFREE(ppiLeafList);
- CKFREE(piLeafCount);
- CKFREE(pdScores);
- FreeProgress(&prProgress);
- CKFREE(pdWeightsL);
- CKFREE(pdWeightsR);
-
-#if TIMING
- StopwatchStop(stopwatch);
- StopwatchDisplay(stdout, "Total time for HHalignWrapper():" , stopwatch);
- StopwatchFree(stopwatch);
-#endif
-
- return dScore; /* FIXME alternative: return averaged pdScores */
-
-}
-/*** end: HHalignWrapper() ***/
git://source.jalview.org / jabaws.git / blobdiff