dist_methods.cpp
上传用户:yhdzpy8989
上传日期:2007-06-13
资源大小:13604k
文件大小:13k
源码类别:

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: dist_methods.cpp,v $
  4.  * PRODUCTION Revision 1000.1  2004/06/01 18:09:24  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.8
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: dist_methods.cpp,v 1000.1 2004/06/01 18:09:24 gouriano Exp $
  10.  * ===========================================================================
  11.  *
  12.  *                            PUBLIC DOMAIN NOTICE
  13.  *               National Center for Biotechnology Information
  14.  *
  15.  *  This software/database is a "United States Government Work" under the
  16.  *  terms of the United States Copyright Act.  It was written as part of
  17.  *  the author's official duties as a United States Government employee and
  18.  *  thus cannot be copyrighted.  This software/database is freely available
  19.  *  to the public for use. The National Library of Medicine and the U.S.
  20.  *  Government have not placed any restriction on its use or reproduction.
  21.  *
  22.  *  Although all reasonable efforts have been taken to ensure the accuracy
  23.  *  and reliability of the software and data, the NLM and the U.S.
  24.  *  Government do not and cannot warrant the performance or results that
  25.  *  may be obtained by using this software or data. The NLM and the U.S.
  26.  *  Government disclaim all warranties, express or implied, including
  27.  *  warranties of performance, merchantability or fitness for any particular
  28.  *  purpose.
  29.  *
  30.  *  Please cite the author in any work or product based on this material.
  31.  *
  32.  * ===========================================================================
  33.  *
  34.  * Authors:  Josh Cherry
  35.  *
  36.  * File Description:  Distance methods for phylogenic tree reconstruction
  37.  *
  38.  */
  39. #include <ncbi_pch.hpp>
  40. #include <algo/phy_tree/dist_methods.hpp>
  41. #include <math.h>
  42. #include "fastme/graph.h"
  43. BEGIN_NCBI_SCOPE
  44. USING_SCOPE(objects);
  45. /// d = -3/4 ln(1 - (4/3)p).
  46. void CDistMethods::JukesCantorDist(const TMatrix& frac_diff,
  47.                                    TMatrix& result)
  48. {
  49.     result.Resize(frac_diff.GetRows(), frac_diff.GetCols());
  50.     for (size_t i = 0;  i < frac_diff.GetRows();  ++i) {
  51.         for (size_t j = 0;  j < frac_diff.GetCols();  ++j) {
  52.             result(i, j) = -0.75 * log(1 - (4.0 / 3.0) * frac_diff(i, j));
  53.         }
  54.     }
  55. }
  56. /// d = -ln(1 - p)
  57. void CDistMethods::PoissonDist(const TMatrix& frac_diff,
  58.                                TMatrix& result)
  59. {
  60.     result.Resize(frac_diff.GetRows(), frac_diff.GetCols());
  61.     for (size_t i = 0;  i < frac_diff.GetRows();  ++i) {
  62.         for (size_t j = 0;  j < frac_diff.GetCols();  ++j) {
  63.             result(i, j) = -log(1 - frac_diff(i, j));
  64.         }
  65.     }
  66. }
  67. /// d = -ln(1 - p - 0.2p^2)
  68. void CDistMethods::KimuraDist(const TMatrix& frac_diff,
  69.                                      TMatrix& result)
  70. {
  71.     result.Resize(frac_diff.GetRows(), frac_diff.GetCols());
  72.     for (size_t i = 0;  i < frac_diff.GetRows();  ++i) {
  73.         for (size_t j = 0;  j < frac_diff.GetCols();  ++j) {
  74.             result(i, j) = -log(1 - frac_diff(i, j)
  75.                                 - 0.2 * frac_diff(i, j) * frac_diff(i, j));
  76.         }
  77.     }
  78. }
  79. /// As per Hillis et al. (Ed.), Molecular Systematics, pg. 488-489
  80. CDistMethods::TTree *CDistMethods::NjTree(const TMatrix& dist_mat,
  81.                                           const vector<string>& labels)
  82. {
  83.     // prepare initial tree (a star phylogeny)
  84.     TTree *tree = new TTree;
  85.     for (unsigned int i = 0;  i < dist_mat.GetRows();  ++i) {
  86.         TTree *new_node = tree->AddNode();
  87.         new_node->GetValue().SetId(i);
  88.         if (labels.empty()) {
  89.             new_node->GetValue().SetLabel() = 'N' + NStr::IntToString(i);
  90.         } else {
  91.             new_node->GetValue().SetLabel() = labels[i];
  92.         }
  93.     }
  94.     // now the real work; do N - 2 neighbor joinings
  95.     int next_id = dist_mat.GetRows();
  96.     TMatrix dmat = dist_mat;
  97.     dmat.Resize(2 * dist_mat.GetRows() - 2, 2 * dist_mat.GetRows() - 2);
  98.     vector<double> r(dmat.GetRows());
  99.     for (unsigned int n = dist_mat.GetRows();  n > 2;  --n) {
  100.         int i, j;
  101.         double m;
  102.         TTree::TNodeList_I neighbor1, neighbor2;
  103.         // first compute r_i
  104.         for (TTree::TNodeList_I it1 = tree->SubNodeBegin();
  105.              it1 != tree->SubNodeEnd();  ++it1) {
  106.             i = (*it1)->GetValue().GetId();
  107.             double sum = 0;
  108.             for (TTree::TNodeList_I it2 = tree->SubNodeBegin();
  109.                  it2 != tree->SubNodeEnd();  ++it2) {
  110.                 if (it1 == it2) {
  111.                     continue;
  112.                 }
  113.                 j = (*it2)->GetValue().GetId();
  114.                 sum += dmat(i, j);
  115.             }
  116.             r[i] = sum;
  117.         }
  118.         // find where M_{i, j} is minimal
  119.         double min_m = 0;
  120.         for (TTree::TNodeList_I it1 = tree->SubNodeBegin();
  121.              it1 != tree->SubNodeEnd();  ++it1) {
  122.             TTree::TNodeList_I it2 = it1;
  123.             ++it2;
  124.             for (;  it2 != tree->SubNodeEnd();  ++it2) {
  125.                 i = (*it1)->GetValue().GetId();
  126.                 j = (*it2)->GetValue().GetId();
  127.                 m = dmat(i, j) - (r[i] + r[j]) / (n - 2);
  128.                 if (m < min_m) {
  129.                     neighbor1 = it1;
  130.                     neighbor2 = it2;
  131.                     min_m = m;
  132.                 }
  133.             }
  134.         }
  135.         // join the neighbors
  136.         TTree *new_node = new TTree;
  137.         i = (*neighbor1)->GetValue().GetId();
  138.         j = (*neighbor2)->GetValue().GetId();
  139.         int u = next_id++;
  140.         new_node->GetValue().SetId(u);
  141.         double viu = dmat(i, j) / 2 + (r[i] - r[j]) / (2 * (n - 2));
  142.         double vju = dmat(i, j) - viu;
  143.         (*neighbor1)->GetValue().SetDist(viu);
  144.         (*neighbor2)->GetValue().SetDist(vju);
  145.         new_node->AddNode(tree->DetachNode(neighbor1));
  146.         new_node->AddNode(tree->DetachNode(neighbor2));
  147.         tree->AddNode(new_node);
  148.         // add new distances to dmat
  149.         for (TTree::TNodeList_CI it1 = tree->SubNodeBegin();
  150.              it1 != tree->SubNodeEnd();  ++it1) {
  151.             int k = (*it1)->GetValue().GetId();
  152.             if (k == u) {
  153.                 continue;
  154.             }
  155.             dmat(k, u) = dmat(u, k)
  156.                 = (dmat(i, k) + dmat(j, k) - dmat(i, j)) / 2;
  157.         }
  158.     }
  159.     // Now the root has just two children, whose distances
  160.     // have not been set.  Could do different things here.
  161.     // Let's make a trifurcation.
  162.     TTree::TNodeList_I it1 = tree->SubNodeBegin();
  163.     TTree::TNodeList_I it2 = it1;
  164.     ++it2;
  165.     if ((*it1)->IsLeaf()) {
  166.         swap(*it1, *it2);
  167.     }
  168.     double d = dmat((*it1)->GetValue().GetId(), (*it2)->GetValue().GetId());
  169.     (*it2)->GetValue().SetDist(d);
  170.     (*it1)->AddNode(tree->DetachNode(it2));
  171.     TTree *rv = tree->DetachNode(it1);
  172.     delete tree;  // just the original root node
  173.     return rv;
  174. }
  175. // The following is a wrapper for Richard Desper's fast minimume evolution
  176. // code.  The user passes in a CDistMethods::TMatrix, and gets back
  177. // a CDistMethods::TTree.  Behind the scenes, the code converts the
  178. // TTMatrix to the representation used by fastme, runs the fastme
  179. // algorithm using fastme_run, and converts the resulting tree
  180. // to the CDistMethods::TTree representation.
  181. BEGIN_SCOPE(fastme)
  182. typedef char boolean;
  183. boolean leaf(meNode *v);
  184. meTree* fastme_run(double** D_in, int numSpecies_in, char **labels,
  185.                    int btype_in, int wtype_in, int ntype_in);
  186. double **initDoubleMatrix(int d);
  187. void freeMatrix(double **D, int size);
  188. void freeTree(meTree *T);
  189. END_SCOPE(fastme)
  190. // Functions to convert from tree representation used by fastme code
  191. static void s_AddFastMeSubtree(fastme::meNode *me_node,
  192.                                CDistMethods::TTree* node,
  193.                                const vector<string>& labels)
  194. {
  195.     if (fastme::leaf(me_node)) {
  196.         int id = NStr::StringToInt(me_node->label);
  197.         node->GetValue().SetId(id);
  198.         if (!labels.empty()) {
  199.             node->GetValue().SetLabel(labels[id]);
  200.         } else {
  201.             node->GetValue().SetLabel(me_node->label);
  202.         }
  203.         return;
  204.     }
  205.     // not a leaf; add two subtrees
  206.     // first left
  207.     CDistMethods::TTree* child_node = node->AddNode();
  208.     child_node->GetValue().SetDist(me_node->leftEdge->distance);
  209.     s_AddFastMeSubtree(me_node->leftEdge->head, child_node, labels);
  210.     // then right
  211.     child_node = node->AddNode();
  212.     child_node->GetValue().SetDist(me_node->rightEdge->distance);
  213.     s_AddFastMeSubtree(me_node->rightEdge->head, child_node, labels);
  214. }
  215. static CDistMethods::TTree* s_ConvertFastMeTree(fastme::meTree *me_tree,
  216.                                                 const vector<string>& labels) {
  217.     CDistMethods::TTree *tree = new CDistMethods::TTree;
  218.     fastme::meEdge *edge;
  219.     edge = me_tree->root->leftEdge;
  220.     CDistMethods::TTree *node = tree->AddNode();
  221.     node->GetValue().SetDist(edge->distance);
  222.     int id = NStr::StringToInt(me_tree->root->label);
  223.     node->GetValue().SetId(id);
  224.     if (!labels.empty()) {
  225.         node->GetValue().SetLabel(labels[id]);
  226.     } else {
  227.         node->GetValue().SetLabel(me_tree->root->label);
  228.     }
  229.     s_AddFastMeSubtree(edge->head, tree, labels);
  230.     return tree;
  231. }
  232. // The publically visible interface
  233. CDistMethods::TTree *CDistMethods::FastMeTree(const TMatrix& dist_mat,
  234.                                               const vector<string>& labels,
  235.                                               EFastMePar btype,
  236.                                               EFastMePar wtype,
  237.                                               EFastMePar ntype)
  238. {
  239.     double **dfme = fastme::initDoubleMatrix(dist_mat.GetRows());
  240.     for (unsigned int i = 0;  i < dist_mat.GetRows();  ++i) {
  241.         for (unsigned int j = 0;  j < dist_mat.GetRows();  ++j) {
  242.             dfme[i][j] = dist_mat(i, j);
  243.         }
  244.     }
  245.     // leaf labels for fastme code; just pass in
  246.     // "0", "1", etc., and fill in real labels later
  247.     char **clabels = new char*[dist_mat.GetRows()];
  248.     vector<string> dummy_labels;
  249.     dummy_labels.resize(dist_mat.GetRows());
  250.     for (unsigned int i = 0;  i < dist_mat.GetRows();  ++i) {
  251.         dummy_labels[i] = NStr::IntToString(i);
  252.         clabels[i] = const_cast<char *>(dummy_labels[i].c_str());
  253.     }
  254.     fastme::meTree *me_out =
  255.         fastme::fastme_run(dfme, dist_mat.GetRows(), clabels,
  256.                            btype, wtype, ntype);
  257.     fastme::freeMatrix(dfme, dist_mat.GetRows());
  258.     delete [] clabels;
  259.     TTree *me_tree = s_ConvertFastMeTree(me_out, labels);
  260.     fastme::freeTree(me_out);
  261.     return me_tree;
  262. }
  263. // Calculate pairwise fractional differences
  264. double CDistMethods::Divergence(const string& seq1, const string& seq2)
  265. {
  266.     int diff_count = 0;
  267.     int pos_count = 0;
  268.     for (unsigned int pos = 0;  pos < seq1.size();  ++pos) {
  269.         if (seq1[pos] == '-' || seq2[pos] == '-') {
  270.             continue;  // ignore this position if a seq has a gap
  271.         }
  272.         pos_count++;
  273.         if (seq1[pos] != seq2[pos]) {
  274.             diff_count++;
  275.         }
  276.     }
  277.     return diff_count / (double) pos_count;
  278. }
  279. void CDistMethods::Divergence(const CAlnVec& avec_in, TMatrix& result)
  280. {
  281.     // want to change gap char of CAlnVec, but no copy constructor,
  282.     // so effectively copy in a round-about way
  283.     CAlnVec avec(avec_in.GetDenseg(), avec_in.GetScope());
  284.     avec.SetGapChar('-');
  285.     int nseqs = avec.GetNumRows();
  286.     result.Resize(nseqs, nseqs);
  287.     string seqi, seqj;
  288.     for (int i = 0;  i < nseqs;  ++i) {
  289.         result(i, i) = 0;  // 0 difference from itself
  290.         avec.GetWholeAlnSeqString(i, seqi);
  291.         for (int j = i + 1;  j < nseqs;  ++j) {
  292.             avec.GetWholeAlnSeqString(j, seqj);
  293.             result(i, j) = result(j, i) = CDistMethods::Divergence(seqi, seqj);
  294.         }
  295.     }
  296. }
  297. #if 0
  298. void CDistMethods::Divergence(const CAlignment& aln, TMatrix& result)
  299. {
  300.     int nseqs = aln.GetSeqs().size();
  301.     result.Resize(nseqs, nseqs);
  302.     for (int i = 0;  i < nseqs;  ++i) {
  303.         result(i, i) = 0;  // 0 difference from itself
  304.         for (int j = i + 1;  j < nseqs;  ++j) {
  305.             result(i, j) = result(j, i) = 
  306.                 CDistMethods::Divergence(aln.GetSeqs()[i], aln.GetSeqs()[j]);
  307.         }
  308.     }
  309. }
  310. #endif
  311. END_NCBI_SCOPE
  312. /*
  313.  * ===========================================================================
  314.  * $Log: dist_methods.cpp,v $
  315.  * Revision 1000.1  2004/06/01 18:09:24  gouriano
  316.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.8
  317.  *
  318.  * Revision 1.8  2004/05/21 21:41:03  gorelenk
  319.  * Added PCH ncbi_pch.hpp
  320.  *
  321.  * Revision 1.7  2004/03/17 17:57:48  jcherry
  322.  * Made fastme set the node ids of leaf nodes
  323.  *
  324.  * Revision 1.6  2004/02/19 16:43:45  jcherry
  325.  * Temporarily disable one form of Divergence() method
  326.  *
  327.  * Revision 1.5  2004/02/19 13:20:04  dicuccio
  328.  * Roll back to version 1.3
  329.  *
  330.  * Revision 1.3  2004/02/10 23:07:46  ucko
  331.  * +<math.h> for log()
  332.  *
  333.  * Revision 1.2  2004/02/10 17:01:19  dicuccio
  334.  * Use swap() instead of iter_swap(), as the latter isn't found on MSVC
  335.  *
  336.  * Revision 1.1  2004/02/10 15:15:58  jcherry
  337.  * Initial version
  338.  *
  339.  * ===========================================================================
  340.  */