binaries/src/clustalo/src/clustal/tree.c

   1 /* -*- mode: c; tab-width: 4; c-basic-offset: 4;  indent-tabs-mode: nil -*- */
   2
   3 /*********************************************************************
   4  * Clustal Omega - Multiple sequence alignment
   5  *
   6  * Copyright (C) 2010 University College Dublin
   7  *
   8  * Clustal-Omega is free software; you can redistribute it and/or
   9  * modify it under the terms of the GNU General Public License as
  10  * published by the Free Software Foundation; either version 2 of the
  11  * License, or (at your option) any later version.
  12  *
  13  * This file is part of Clustal-Omega.
  14  *
  15  ********************************************************************/
  16
  17 /*
  18  *  RCS $Id: tree.c 278 2013-05-16 15:53:45Z fabian $
  19  */
  20
  21 #ifdef HAVE_CONFIG_H
  22 #include "config.h"
  23 #endif
  24
  25 #include <stdlib.h>
  26 #include <string.h>
  27 #include <assert.h>
  28
  29 #include "util.h"
  30 #include "log.h"
  31 #include "muscle_upgma.h"
  32 #include "tree.h"
  33
  34 /**
  35  *
  36  * @brief Creates a UPGMA guide tree. This is a frontend function to
  37  * the ported Muscle UPGMA code ().
  38  *
  39  * @param[out] tree
  40  * created upgma tree. will be allocated here. use FreeMuscleTree()
  41  * to free
  42  * @param[in] labels
  43  * pointer to nseq sequence names
  44  * @param[in] distmat
  45  * distance matrix
  46  * @param[in] ftree
  47  * optional: if non-NULL, tree will be written to this files
  48  *
  49  * @see FreeMuscleTree()
  50  * @see MuscleUpgma2()
  51  *
  52  */
  53 void
  54 GuideTreeUpgma(tree_t **tree, char **labels,
  55                 symmatrix_t *distmat, char *ftree)
  56 {
  57     linkage_t linkage = LINKAGE_AVG;
  58     FILE *fp = NULL;
  59
  60     if (NULL != ftree) {
  61         if (NULL == (fp=fopen(ftree, "w"))) {
  62             Log(&rLog, LOG_ERROR, "Couldn't open tree-file '%s' for writing. Skipping", ftree);
  63         }
  64         /* fp NULL is handled later */
  65     }
  66
  67     (*tree) = (tree_t *) CKMALLOC(1 * sizeof(tree_t));
  68     MuscleUpgma2((*tree), distmat, linkage, labels);
  69
  70     if (rLog.iLogLevelEnabled <= LOG_DEBUG) {
  71         Log(&rLog, LOG_DEBUG, "tree logging...");
  72         LogTree((*tree), LogGetFP(&rLog, LOG_DEBUG));
  73     }
  74
  75     if (NULL != fp) {
  76         MuscleTreeToFile(fp, (*tree));
  77         Log(&rLog, LOG_INFO, "Guide tree written to %s", ftree);
  78         fclose(fp);
  79     }
  80 }
  81 /***   end: guidetree_upgma   ***/
  82
  83
  84
  85 /**
  86  *
  87  * @brief
  88  *
  89  * @param[out] tree
  90  * created upgma tree. will be allocated here. use FreeMuscleTree()
  91  * to free
  92  * @param[in] mseq
  93  * @param[in] ftree
  94  *
  95  * @return non-zero on error
  96  *
  97  */
  98 int
  99 GuideTreeFromFile(tree_t **tree, mseq_t *mseq, char *ftree)
 100 {
 101     int iNodeCount;
 102     int iNodeIndex;
 103
 104     (*tree) = (tree_t *) CKMALLOC(1 * sizeof(tree_t));
 105     if (MuscleTreeFromFile((*tree), ftree)!=0) {
 106         Log(&rLog, LOG_ERROR, "%s", "MuscleTreeFromFile failed");
 107         return -1;
 108     }
 109
 110     /* Make sure tree is rooted */
 111     if (!IsRooted((*tree))) {
 112         Log(&rLog, LOG_ERROR, "User tree must be rooted");
 113         return -1;
 114     }
 115
 116     if ((int)GetLeafCount((*tree)) != mseq->nseqs) {
 117         Log(&rLog, LOG_ERROR, "User tree does not match input sequences");
 118         return -1;
 119     }
 120
 121     /* compare tree labels and sequence names and set leaf-ids */
 122     iNodeCount = GetNodeCount((*tree));
 123     for (iNodeIndex = 0; iNodeIndex < iNodeCount; ++iNodeIndex) {
 124         char *LeafName;
 125         int iSeqIndex;
 126
 127         if (!IsLeaf(iNodeIndex, (*tree)))
 128             continue;
 129         LeafName = GetLeafName(iNodeIndex, (*tree));
 130
 131         if ((iSeqIndex=FindSeqName(LeafName, mseq))==-1) {
 132             Log(&rLog, LOG_ERROR, "Label '%s' in tree could not be found in sequence names", LeafName);
 133             return -1;
 134         }
 135
 136         SetLeafId((*tree), iNodeIndex, iSeqIndex);
 137     }
 138
 139     if (rLog.iLogLevelEnabled <= LOG_DEBUG) {
 140         Log(&rLog, LOG_DEBUG, "tree logging...");
 141         LogTree((*tree),  LogGetFP(&rLog, LOG_DEBUG));
 142     }
 143
 144     return 0;
 145 }
 146 /***   end: GuideTreeFromFile()   ***/
 147
 148
 149
 150 /**
 151  *
 152  * @brief Depth first traversal of tree, i.e. leaf nodes (sequences)
 153  * will be visited first. Order can be used to guide progressive
 154  * alignment order.
 155  *
 156  * @param[out] piOrderLR_p
 157  * order in which left/right nodes (profiles) are to be aligned.
 158  * allocated here; caller must free.
 159  * @param[in] tree
 160  * The tree to traverse; has to be rooted
 161  * @param[in] mseq
 162  * corresponding multiple sequence structure
 163  *
 164  */
 165 void
 166 TraverseTree(int **piOrderLR_p,
 167               tree_t *tree, mseq_t *mseq)
 168 {
 169     int tree_nodeindex = 0;
 170     int order_index = 0;
 171     int iLeafCount = 0;
 172
 173     assert(NULL!=tree);
 174     assert(NULL!=mseq);
 175     assert(IsRooted(tree));
 176
 177     /* allocate memory for node/profile alignment order;
 178      * for every node allocate DIFF_NODE (3) int (1 left, 1 right, 1 parent)
 179      */
 180     *piOrderLR_p = (int *)CKCALLOC(DIFF_NODE * GetNodeCount(tree), sizeof(int));
 181
 182     /* Log(&rLog, LOG_FORCED_DEBUG, "print tree->m_iNodeCount=%d", tree->m_iNodeCount); */
 183
 184
 185     tree_nodeindex = FirstDepthFirstNode(tree);
 186     /*LOG_DEBUG("Starting with treenodeindex = %d", tree_nodeindex);*/
 187
 188     order_index = 0;
 189
 190     do {
 191         if (IsLeaf(tree_nodeindex, tree)) {
 192             int leafid = GetLeafId(tree_nodeindex, tree);
 193             if (leafid >= mseq->nseqs){
 194                 Log(&rLog, LOG_FATAL, "Sequence index out of range during tree traversal (leafid=%d nseqs=%d)",
 195                     leafid, mseq->nseqs);
 196             }
 197             if (NULL != mseq->tree_order){
 198                 mseq->tree_order[iLeafCount] = leafid;
 199                 iLeafCount++;
 200             }
 201
 202             /* this is a leaf node,
 203              * indicate this by registering same leafid for left/right
 204              */
 205
 206             (*piOrderLR_p)[DIFF_NODE*order_index+LEFT_NODE] = leafid;
 207             (*piOrderLR_p)[DIFF_NODE*order_index+RGHT_NODE] = leafid;
 208             (*piOrderLR_p)[DIFF_NODE*order_index+PRNT_NODE] = tree_nodeindex;
 209
 210             Log(&rLog, LOG_DEBUG, "Tree traversal: Visited leaf-node %d (leaf-id %d = Seq '%s')",
 211                  tree_nodeindex, leafid, mseq->sqinfo[leafid].name);
 212
 213         } else {
 214             int merge_nodeindex;
 215             int left;
 216             int right;
 217
 218             merge_nodeindex = tree_nodeindex;
 219             left  = GetLeft(tree_nodeindex, tree);
 220             right = GetRight(tree_nodeindex, tree);
 221
 222             /* this is not a leaf node but a merge node,
 223              * register left node (even) and right node (odd)
 224              */
 225             (*piOrderLR_p)[DIFF_NODE*order_index+LEFT_NODE] = left;
 226             (*piOrderLR_p)[DIFF_NODE*order_index+RGHT_NODE] = right;
 227             (*piOrderLR_p)[DIFF_NODE*order_index+PRNT_NODE] = merge_nodeindex;
 228
 229             Log(&rLog, LOG_DEBUG, "Tree traversal: Visited non-leaf node %d with siblings %d (L) and %d (R)",
 230                  merge_nodeindex, left, right);
 231         }
 232         tree_nodeindex = NextDepthFirstNode(tree_nodeindex, tree);
 233
 234         order_index++;
 235
 236     } while (NULL_NEIGHBOR != tree_nodeindex);
 237
 238     return;
 239 }
 240 /***   end: TraverseTree   ***/