+++ /dev/null
-/////////////////////////////////////////////////////////////////
-// EvolutionaryTree.hpp
-//
-// Utilities for reading/writing multiple sequence data.
-/////////////////////////////////////////////////////////////////
-
-#ifndef __EVOLUTIONARYTREE_HPP__
-#define __EVOLUTIONARYTREE_HPP__
-
-#include <string>
-#include <list>
-#include <stdio.h>
-#include "SafeVector.h"
-#include "MultiSequence.h"
-#include "Sequence.h"
-#include "Util.hpp"
-
-using namespace std;
-
-
-/////////////////////////////////////////////////////////////////
-// TreeNode
-//
-// The fundamental unit for representing an alignment tree. The
-// guide tree is represented as a binary tree.
-/////////////////////////////////////////////////////////////////
-namespace MXSCARNA {
-class TreeNode {
- int sequenceLabel; // sequence label
- float sequenceIdentity; // sequence identity
- TreeNode *left, *right, *parent; // pointers to left, right children
- float leftLength, rightLength; // the length of left and right edge
- /////////////////////////////////////////////////////////////////
- // TreeNode::PrintNode()
- //
- // Internal routine used to print out the sequence comments
- // associated with the evolutionary tree, using a hierarchical
- // parenthesized format.
- /////////////////////////////////////////////////////////////////
-
- void PrintNode (ostream &outfile, const MultiSequence *sequences) const {
-
- // if this is a leaf node, print out the associated sequence comment
- if (sequenceLabel >= 0)
- //outfile << sequences->GetSequence (sequenceLabel)->GetHeader();
- outfile << sequences->GetSequence (sequenceLabel)->GetLabel();
-
- // otherwise, it must have two children; print out their subtrees recursively
- else {
- assert (left);
- assert (right);
-
- outfile << "(";
- left->PrintNode (outfile, sequences);
- outfile << ",";
- right->PrintNode (outfile, sequences);
- outfile << ")";
- }
- }
-
- public:
-
- /////////////////////////////////////////////////////////////////
- // TreeNode::TreeNode()
- //
- // Constructor for a tree node. Note that sequenceLabel = -1
- // implies that the current node is not a leaf in the tree.
- /////////////////////////////////////////////////////////////////
-
- TreeNode (int sequenceLabel) : sequenceLabel (sequenceLabel),
- left (NULL), right (NULL), parent (NULL) {
- assert (sequenceLabel >= -1);
- }
-
- /////////////////////////////////////////////////////////////////
- // TreeNode::~TreeNode()
- //
- // Destructor for a tree node. Recursively deletes all children.
- /////////////////////////////////////////////////////////////////
-
- ~TreeNode (){
- if (left){ delete left; left = NULL; }
- if (right){ delete right; right = NULL; }
- parent = NULL;
- }
-
-
- // getters
- int GetSequenceLabel () const { return sequenceLabel; }
- TreeNode *GetLeftChild () const { return left; }
- TreeNode *GetRightChild () const { return right; }
- TreeNode *GetParent () const { return parent; }
- float GetIdentity () const { return sequenceIdentity; }
- float GetLeftLength () const { return leftLength; }
- float GetRightLength () const { return rightLength; }
- // setters
- void SetSequenceLabel (int sequenceLabel){ this->sequenceLabel = sequenceLabel; assert (sequenceLabel >= -1); }
- void SetLeftChild (TreeNode *left){ this->left = left; }
- void SetRightChild (TreeNode *right){ this->right = right; }
- void SetParent (TreeNode *parent){ this->parent = parent; }
- void SetIdentity (float identity) { this->sequenceIdentity = identity; }
- void SetLeftLength (float identity) { this->leftLength = identity; }
- void SetRightLength (float identity) {this->rightLength = identity; }
- /////////////////////////////////////////////////////////////////
- // TreeNode::ComputeTree()
- //
- // Routine used to compute an evolutionary tree based on the
- // given distance matrix. We assume the distance matrix has the
- // form, distMatrix[i][j] = expected accuracy of aligning i with j.
- /////////////////////////////////////////////////////////////////
-
- static TreeNode *ComputeTree (const VVF &distMatrix, const VVF &identityMatrix){
-
- int numSeqs = distMatrix.size(); // number of sequences in distance matrix
- VVF distances (numSeqs, VF (numSeqs)); // a copy of the distance matrix
- SafeVector<TreeNode *> nodes (numSeqs, NULL); // list of nodes for each sequence
- SafeVector<int> valid (numSeqs, 1); // valid[i] tells whether or not the ith
- // nodes in the distances and nodes array
- // are valid
- VVF identities (numSeqs, VF (numSeqs));
- SafeVector<int> countCluster (numSeqs, 1);
-
- // initialization: make a copy of the distance matrix
- for (int i = 0; i < numSeqs; i++) {
- for (int j = 0; j < numSeqs; j++) {
- distances[i][j] = distMatrix[i][j];
- identities[i][j] = identityMatrix[i][j];
- }
- }
-
- // initialization: create all the leaf nodes
- for (int i = 0; i < numSeqs; i++){
- nodes[i] = new TreeNode (i);
- assert (nodes[i]);
- }
-
- // repeat until only a single node left
- for (int numNodesLeft = numSeqs; numNodesLeft > 1; numNodesLeft--){
- float bestProb = -1;
- pair<int,int> bestPair;
-
- // find the closest pair
- for (int i = 0; i < numSeqs; i++) if (valid[i]){
- for (int j = i+1; j < numSeqs; j++) if (valid[j]){
- if (distances[i][j] > bestProb){
- bestProb = distances[i][j];
- bestPair = make_pair(i, j);
- }
- }
- }
-
- // merge the closest pair
- TreeNode *newParent = new TreeNode (-1);
- newParent->SetLeftChild (nodes[bestPair.first]);
- newParent->SetRightChild (nodes[bestPair.second]);
- nodes[bestPair.first]->SetParent (newParent);
- nodes[bestPair.second]->SetParent (newParent);
- nodes[bestPair.first] = newParent;
- nodes[bestPair.second] = NULL;
- newParent->SetIdentity(identities[bestPair.first][bestPair.second]);
-
-
- // now update the distance matrix
- for (int i = 0; i < numSeqs; i++) if (valid[i]){
- distances[bestPair.first][i] = distances[i][bestPair.first]
- = (distances[i][bestPair.first]*countCluster[bestPair.first]
- + distances[i][bestPair.second]*countCluster[bestPair.second])
- / (countCluster[bestPair.first] + countCluster[bestPair.second]);
-// distances[bestPair.first][i] = distances[i][bestPair.first]
-// = (distances[i][bestPair.first] + distances[i][bestPair.second]) * bestProb / 2;
- identities[bestPair.first][i] = identities[i][bestPair.first]
- = (identities[i][bestPair.first]*countCluster[bestPair.first]
- + identities[i][bestPair.second]*countCluster[bestPair.second])
- / (countCluster[bestPair.first] + countCluster[bestPair.second]);
- }
-
- // finally, mark the second node entry as no longer valid
- countCluster[bestPair.first] += countCluster[bestPair.second];
- valid[bestPair.second] = 0;
- }
-
- assert (nodes[0]);
- return nodes[0];
- }
-
- /////////////////////////////////////////////////////////////////
- // TreeNode::Print()
- //
- // Print out the subtree associated with this node in a
- // parenthesized representation.
- /////////////////////////////////////////////////////////////////
-
- void Print (ostream &outfile, const MultiSequence *sequences) const {
-// outfile << "Alignment tree: ";
- PrintNode (outfile, sequences);
- outfile << endl;
- }
-};
-}
-#endif //__EVOLUTIONARYTREE_HPP__
git://source.jalview.org / jabaws.git / blobdiff