website/archive/binaries/mac/src/muscle/progalign.cpp

   1 #include "muscle.h"\r
   2 #include "tree.h"\r
   3 #include "seqvect.h"\r
   4 #include "profile.h"\r
   5 #include "msa.h"\r
   6 #include "pwpath.h"\r
   7 #include "distfunc.h"\r
   8 #include "textfile.h"\r
   9 #include "estring.h"\r
  10 \r
  11 #define TRACE           0\r
  12 #define VALIDATE        0\r
  13 #define TRACE_LENGTH_DELTA      0\r
  14 \r
  15 static void LogLeafNames(const Tree &tree, unsigned uNodeIndex)\r
  16         {\r
  17         const unsigned uNodeCount = tree.GetNodeCount();\r
  18         unsigned *Leaves = new unsigned[uNodeCount];\r
  19         unsigned uLeafCount;\r
  20         GetLeaves(tree, uNodeIndex, Leaves, &uLeafCount);\r
  21         for (unsigned i = 0; i < uLeafCount; ++i)\r
  22                 {\r
  23                 if (i > 0)\r
  24                         Log(",");\r
  25                 Log("%s", tree.GetLeafName(Leaves[i]));\r
  26                 }\r
  27         delete[] Leaves;\r
  28         }\r
  29 \r
  30 ProgNode *ProgressiveAlignE(const SeqVect &v, const Tree &GuideTree, MSA &a)\r
  31         {\r
  32         assert(GuideTree.IsRooted());\r
  33 \r
  34 #if     TRACE\r
  35         Log("GuideTree:\n");\r
  36         GuideTree.LogMe();\r
  37 #endif\r
  38 \r
  39         const unsigned uSeqCount = v.Length();\r
  40         const unsigned uNodeCount = 2*uSeqCount - 1;\r
  41         const unsigned uIterCount = uSeqCount - 1;\r
  42 \r
  43         WEIGHT *Weights = new WEIGHT[uSeqCount];\r
  44         CalcClustalWWeights(GuideTree, Weights);\r
  45 \r
  46         ProgNode *ProgNodes = new ProgNode[uNodeCount];\r
  47 \r
  48         unsigned uJoin = 0;\r
  49         unsigned uTreeNodeIndex = GuideTree.FirstDepthFirstNode();\r
  50         SetProgressDesc("Align node");\r
  51         do\r
  52                 {\r
  53                 if (GuideTree.IsLeaf(uTreeNodeIndex))\r
  54                         {\r
  55                         if (uTreeNodeIndex >= uNodeCount)\r
  56                                 Quit("TreeNodeIndex=%u NodeCount=%u\n", uTreeNodeIndex, uNodeCount);\r
  57                         ProgNode &Node = ProgNodes[uTreeNodeIndex];\r
  58                         unsigned uId = GuideTree.GetLeafId(uTreeNodeIndex);\r
  59                         if (uId >= uSeqCount)\r
  60                                 Quit("Seq index out of range");\r
  61                         const Seq &s = *(v[uId]);\r
  62                         Node.m_MSA.FromSeq(s);\r
  63                         Node.m_MSA.SetSeqId(0, uId);\r
  64                         Node.m_uLength = Node.m_MSA.GetColCount();\r
  65                         Node.m_Weight = Weights[uId];\r
  66                 // TODO: Term gaps settable\r
  67                         Node.m_Prof = ProfileFromMSA(Node.m_MSA);\r
  68                         Node.m_EstringL = 0;\r
  69                         Node.m_EstringR = 0;\r
  70 #if     TRACE\r
  71                         Log("Leaf id=%u\n", uId);\r
  72                         Log("MSA=\n");\r
  73                         Node.m_MSA.LogMe();\r
  74                         Log("Profile (from MSA)=\n");\r
  75                         ListProfile(Node.m_Prof, Node.m_uLength, &Node.m_MSA);\r
  76 #endif\r
  77                         }\r
  78                 else\r
  79                         {\r
  80                         Progress(uJoin, uSeqCount - 1);\r
  81                         ++uJoin;\r
  82 \r
  83                         const unsigned uMergeNodeIndex = uTreeNodeIndex;\r
  84                         ProgNode &Parent = ProgNodes[uMergeNodeIndex];\r
  85 \r
  86                         const unsigned uLeft = GuideTree.GetLeft(uTreeNodeIndex);\r
  87                         const unsigned uRight = GuideTree.GetRight(uTreeNodeIndex);\r
  88 \r
  89                         if (g_bVerbose)\r
  90                                 {\r
  91                                 Log("Align: (");\r
  92                                 LogLeafNames(GuideTree, uLeft);\r
  93                                 Log(") (");\r
  94                                 LogLeafNames(GuideTree, uRight);\r
  95                                 Log(")\n");\r
  96                                 }\r
  97 \r
  98                         ProgNode &Node1 = ProgNodes[uLeft];\r
  99                         ProgNode &Node2 = ProgNodes[uRight];\r
 100 \r
 101 #if     TRACE\r
 102                         Log("AlignTwoMSAs:\n");\r
 103 #endif\r
 104                         AlignTwoProfs(\r
 105                           Node1.m_Prof, Node1.m_uLength, Node1.m_Weight,\r
 106                           Node2.m_Prof, Node2.m_uLength, Node2.m_Weight,\r
 107                           Parent.m_Path,\r
 108                           &Parent.m_Prof, &Parent.m_uLength);\r
 109 #if     TRACE_LENGTH_DELTA\r
 110                         {\r
 111                         unsigned L = Node1.m_uLength;\r
 112                         unsigned R = Node2.m_uLength;\r
 113                         unsigned P = Parent.m_Path.GetEdgeCount();\r
 114                         unsigned Max = L > R ? L : R;\r
 115                         unsigned d = P - Max;\r
 116                         Log("LD%u;%u;%u;%u\n", L, R, P, d);\r
 117                         }\r
 118 #endif\r
 119                         PathToEstrings(Parent.m_Path, &Parent.m_EstringL, &Parent.m_EstringR);\r
 120 \r
 121                         Parent.m_Weight = Node1.m_Weight + Node2.m_Weight;\r
 122 \r
 123 #if     VALIDATE\r
 124                         {\r
 125 #if     TRACE\r
 126                         Log("AlignTwoMSAs:\n");\r
 127 #endif\r
 128                         PWPath TmpPath;\r
 129                         AlignTwoMSAs(Node1.m_MSA, Node2.m_MSA, Parent.m_MSA, TmpPath);\r
 130                         ProfPos *P1 = ProfileFromMSA(Node1.m_MSA, true);\r
 131                         ProfPos *P2 = ProfileFromMSA(Node2.m_MSA, true);\r
 132                         unsigned uLength = Parent.m_MSA.GetColCount();\r
 133                         ProfPos *TmpProf = ProfileFromMSA(Parent.m_MSA, true);\r
 134 \r
 135 #if     TRACE\r
 136                         Log("Node1 MSA=\n");\r
 137                         Node1.m_MSA.LogMe();\r
 138 \r
 139                         Log("Node1 prof=\n");\r
 140                         ListProfile(Node1.m_Prof, Node1.m_MSA.GetColCount(), &Node1.m_MSA);\r
 141                         Log("Node1 prof (from MSA)=\n");\r
 142                         ListProfile(P1, Node1.m_MSA.GetColCount(), &Node1.m_MSA);\r
 143 \r
 144                         AssertProfsEq(Node1.m_Prof, Node1.m_uLength, P1, Node1.m_MSA.GetColCount());\r
 145 \r
 146                         Log("Node2 prof=\n");\r
 147                         ListProfile(Node2.m_Prof, Node2.m_MSA.GetColCount(), &Node2.m_MSA);\r
 148 \r
 149                         Log("Node2 MSA=\n");\r
 150                         Node2.m_MSA.LogMe();\r
 151 \r
 152                         Log("Node2 prof (from MSA)=\n");\r
 153                         ListProfile(P2, Node2.m_MSA.GetColCount(), &Node2.m_MSA);\r
 154 \r
 155                         AssertProfsEq(Node2.m_Prof, Node2.m_uLength, P2, Node2.m_MSA.GetColCount());\r
 156 \r
 157                         TmpPath.AssertEqual(Parent.m_Path);\r
 158 \r
 159                         Log("Parent MSA=\n");\r
 160                         Parent.m_MSA.LogMe();\r
 161 \r
 162                         Log("Parent prof=\n");\r
 163                         ListProfile(Parent.m_Prof, Parent.m_uLength, &Parent.m_MSA);\r
 164 \r
 165                         Log("Parent prof (from MSA)=\n");\r
 166                         ListProfile(TmpProf, Parent.m_MSA.GetColCount(), &Parent.m_MSA);\r
 167 \r
 168 #endif  // TRACE\r
 169                         AssertProfsEq(Parent.m_Prof, Parent.m_uLength,\r
 170                           TmpProf, Parent.m_MSA.GetColCount());\r
 171                         delete[] P1;\r
 172                         delete[] P2;\r
 173                         delete[] TmpProf;\r
 174                         }\r
 175 #endif  // VALIDATE\r
 176 \r
 177                         Node1.m_MSA.Clear();\r
 178                         Node2.m_MSA.Clear();\r
 179 \r
 180                 // Don't delete profiles, may need them for tree refinement.\r
 181                         //delete[] Node1.m_Prof;\r
 182                         //delete[] Node2.m_Prof;\r
 183                         //Node1.m_Prof = 0;\r
 184                         //Node2.m_Prof = 0;\r
 185                         }\r
 186                 uTreeNodeIndex = GuideTree.NextDepthFirstNode(uTreeNodeIndex);\r
 187                 }\r
 188         while (NULL_NEIGHBOR != uTreeNodeIndex);\r
 189         ProgressStepsDone();\r
 190 \r
 191         if (g_bBrenner)\r
 192                 MakeRootMSABrenner((SeqVect &) v, GuideTree, ProgNodes, a);\r
 193         else\r
 194                 MakeRootMSA(v, GuideTree, ProgNodes, a);\r
 195 \r
 196 #if     VALIDATE\r
 197         {\r
 198         unsigned uRootNodeIndex = GuideTree.GetRootNodeIndex();\r
 199         const ProgNode &RootProgNode = ProgNodes[uRootNodeIndex];\r
 200         AssertMSAEq(a, RootProgNode.m_MSA);\r
 201         }\r
 202 #endif\r
 203 \r
 204         delete[] Weights;\r
 205         return ProgNodes;\r
 206         }\r