binaries/src/muscle/muscle.cpp

   1 #include "muscle.h"\r
   2 #include "msa.h"\r
   3 #include "seqvect.h"\r
   4 #include "msa.h"\r
   5 #include "tree.h"\r
   6 #include "profile.h"\r
   7 \r
   8 void MUSCLE(SeqVect &v, MSA &msaOut)\r
   9         {\r
  10         const unsigned uSeqCount = v.Length();\r
  11 \r
  12         if (0 == uSeqCount)\r
  13                 Quit("No sequences in input file");\r
  14 \r
  15         ALPHA Alpha = ALPHA_Undefined;\r
  16         switch (g_SeqType)\r
  17                 {\r
  18         case SEQTYPE_Auto:\r
  19                 Alpha = v.GuessAlpha();\r
  20                 break;\r
  21 \r
  22         case SEQTYPE_Protein:\r
  23                 Alpha = ALPHA_Amino;\r
  24                 break;\r
  25 \r
  26         case SEQTYPE_RNA:\r
  27                 Alpha = ALPHA_RNA;\r
  28                 break;\r
  29 \r
  30         case SEQTYPE_DNA:\r
  31                 Alpha = ALPHA_DNA;\r
  32                 break;\r
  33 \r
  34         default:\r
  35                 Quit("Invalid seq type");\r
  36                 }\r
  37         SetAlpha(Alpha);\r
  38         v.FixAlpha();\r
  39 \r
  40         if (ALPHA_DNA == Alpha || ALPHA_RNA == Alpha)\r
  41                 {\r
  42                 SetPPScore(PPSCORE_SPN);\r
  43                 g_Distance1 = DISTANCE_Kmer4_6;\r
  44                 }\r
  45 \r
  46         unsigned uMaxL = 0;\r
  47         unsigned uTotL = 0;\r
  48         for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
  49                 {\r
  50                 unsigned L = v.GetSeq(uSeqIndex).Length();\r
  51                 uTotL += L;\r
  52                 if (L > uMaxL)\r
  53                         uMaxL = L;\r
  54                 }\r
  55 \r
  56         SetIter(1);\r
  57         g_bDiags = g_bDiags1;\r
  58         SetSeqStats(uSeqCount, uMaxL, uTotL/uSeqCount);\r
  59 \r
  60         MSA::SetIdCount(uSeqCount);\r
  61 \r
  62 //// Initialize sequence ids.\r
  63 //// From this point on, ids must somehow propogate from here.\r
  64 //      for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
  65 //              v.SetSeqId(uSeqIndex, uSeqIndex);\r
  66 \r
  67         if (uSeqCount > 1)\r
  68                 MHackStart(v);\r
  69 \r
  70         if (0 == uSeqCount)\r
  71                 {\r
  72                 msaOut.Clear();\r
  73                 return;\r
  74                 }\r
  75 \r
  76         if (1 == uSeqCount && ALPHA_Amino == Alpha)\r
  77                 {\r
  78                 const Seq &s = v.GetSeq(0);\r
  79                 msaOut.FromSeq(s);\r
  80                 return;\r
  81                 }\r
  82 \r
  83 // First iteration\r
  84         Tree GuideTree;\r
  85         TreeFromSeqVect(v, GuideTree, g_Cluster1, g_Distance1, g_Root1);\r
  86 \r
  87         SetMuscleTree(GuideTree);\r
  88 \r
  89         ProgNode *ProgNodes = 0;\r
  90         if (g_bLow)\r
  91                 ProgNodes = ProgressiveAlignE(v, GuideTree, msaOut);\r
  92         else\r
  93                 ProgressiveAlign(v, GuideTree, msaOut);\r
  94         SetCurrentAlignment(msaOut);\r
  95 \r
  96         if (1 == g_uMaxIters || 2 == uSeqCount)\r
  97                 {\r
  98                 MHackEnd(msaOut);\r
  99                 return;\r
 100                 }\r
 101 \r
 102         g_bDiags = g_bDiags2;\r
 103         SetIter(2);\r
 104 \r
 105         if (g_bLow)\r
 106                 {\r
 107                 if (0 != g_uMaxTreeRefineIters)\r
 108                         RefineTreeE(msaOut, v, GuideTree, ProgNodes);\r
 109                 }\r
 110         else\r
 111                 RefineTree(msaOut, GuideTree);\r
 112 \r
 113         extern void DeleteProgNode(ProgNode &Node);\r
 114         const unsigned uNodeCount = GuideTree.GetNodeCount();\r
 115         for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)\r
 116                 DeleteProgNode(ProgNodes[uNodeIndex]);\r
 117 \r
 118         delete[] ProgNodes;\r
 119         ProgNodes = 0;\r
 120 \r
 121         SetSeqWeightMethod(g_SeqWeight2);\r
 122         SetMuscleTree(GuideTree);\r
 123 \r
 124         if (g_bAnchors)\r
 125                 RefineVert(msaOut, GuideTree, g_uMaxIters - 2);\r
 126         else\r
 127                 RefineHoriz(msaOut, GuideTree, g_uMaxIters - 2, false, false);\r
 128 \r
 129         MHackEnd(msaOut);\r
 130         }\r