binaries/src/muscle/profile.cpp

   1 #include "muscle.h"\r
   2 #include "textfile.h"\r
   3 #include "msa.h"\r
   4 #include "tree.h"\r
   5 #include "profile.h"\r
   6 #include "objscore.h"\r
   7 \r
   8 bool TreeNeededForWeighting(SEQWEIGHT s)\r
   9         {\r
  10         switch (s)\r
  11                 {\r
  12         case SEQWEIGHT_ClustalW:\r
  13         case SEQWEIGHT_ThreeWay:\r
  14                 return true;\r
  15         default:\r
  16                 return false;\r
  17                 }\r
  18         }\r
  19 \r
  20 static ProfPos *ProfileFromMSALocal(MSA &msa, Tree &tree)\r
  21         {\r
  22         const unsigned uSeqCount = msa.GetSeqCount();\r
  23         for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
  24                 msa.SetSeqId(uSeqIndex, uSeqIndex);\r
  25 \r
  26         if (TreeNeededForWeighting(g_SeqWeight2))\r
  27                 {\r
  28                 TreeFromMSA(msa, tree, g_Cluster2, g_Distance2, g_Root1);\r
  29                 SetMuscleTree(tree);\r
  30                 }\r
  31         return ProfileFromMSA(msa);\r
  32         }\r
  33 \r
  34 void ProfileProfile(MSA &msa1, MSA &msa2, MSA &msaOut)\r
  35         {\r
  36         //ALPHA Alpha = ALPHA_Undefined;\r
  37         //switch (g_SeqType)\r
  38         //      {\r
  39         //case SEQTYPE_Auto:\r
  40         //      Alpha = msa1.GuessAlpha();\r
  41         //      break;\r
  42 \r
  43         //case SEQTYPE_Protein:\r
  44         //      Alpha = ALPHA_Amino;\r
  45         //      break;\r
  46 \r
  47         //case SEQTYPE_DNA:\r
  48         //      Alpha = ALPHA_DNA;\r
  49         //      break;\r
  50 \r
  51         //case SEQTYPE_RNA:\r
  52         //      Alpha = ALPHA_RNA;\r
  53         //      break;\r
  54 \r
  55         //default:\r
  56         //      Quit("Invalid SeqType");\r
  57         //      }\r
  58         //SetAlpha(Alpha);\r
  59 \r
  60         //msa1.FixAlpha();\r
  61         //msa2.FixAlpha();\r
  62 \r
  63         unsigned uLength1;\r
  64         unsigned uLength2;\r
  65 \r
  66         uLength1 = msa1.GetColCount();\r
  67         uLength2 = msa2.GetColCount();\r
  68 \r
  69         Tree tree1;\r
  70         Tree tree2;\r
  71         ProfPos *Prof1 = ProfileFromMSALocal(msa1, tree1);\r
  72         ProfPos *Prof2 = ProfileFromMSALocal(msa2, tree2);\r
  73 \r
  74         PWPath Path;\r
  75         ProfPos *ProfOut;\r
  76         unsigned uLengthOut;\r
  77         Progress("Aligning profiles");\r
  78         AlignTwoProfs(Prof1, uLength1, 1.0, Prof2, uLength2, 1.0, Path, &ProfOut, &uLengthOut);\r
  79 \r
  80         Progress("Building output");\r
  81         AlignTwoMSAsGivenPath(Path, msa1, msa2, msaOut);\r
  82         }\r
  83 \r
  84 // Do profile-profile alignment\r
  85 void Profile()\r
  86         {\r
  87         if (0 == g_pstrFileName1 || 0 == g_pstrFileName2)\r
  88                 Quit("-profile needs -in1 and -in2");\r
  89 \r
  90         SetSeqWeightMethod(g_SeqWeight1);\r
  91 \r
  92         TextFile file1(g_pstrFileName1);\r
  93         TextFile file2(g_pstrFileName2);\r
  94 \r
  95         MSA msa1;\r
  96         MSA msa2;\r
  97         MSA msaOut;\r
  98 \r
  99         Progress("Reading %s", g_pstrFileName1);\r
 100         msa1.FromFile(file1);\r
 101         Progress("%u seqs %u cols", msa1.GetSeqCount(), msa1.GetColCount());\r
 102 \r
 103         Progress("Reading %s", g_pstrFileName2);\r
 104         msa2.FromFile(file2);\r
 105         Progress("%u seqs %u cols", msa2.GetSeqCount(), msa2.GetColCount());\r
 106 \r
 107         ALPHA Alpha = ALPHA_Undefined;\r
 108         switch (g_SeqType)\r
 109                 {\r
 110         case SEQTYPE_Auto:\r
 111                 Alpha = msa1.GuessAlpha();\r
 112                 break;\r
 113 \r
 114         case SEQTYPE_Protein:\r
 115                 Alpha = ALPHA_Amino;\r
 116                 break;\r
 117 \r
 118         case SEQTYPE_DNA:\r
 119                 Alpha = ALPHA_DNA;\r
 120                 break;\r
 121 \r
 122         case SEQTYPE_RNA:\r
 123                 Alpha = ALPHA_RNA;\r
 124                 break;\r
 125 \r
 126         default:\r
 127                 Quit("Invalid seq type");\r
 128                 }\r
 129         SetAlpha(Alpha);\r
 130 \r
 131         msa1.FixAlpha();\r
 132         msa2.FixAlpha();\r
 133 \r
 134         SetPPScore();\r
 135         if (ALPHA_DNA == Alpha || ALPHA_RNA == Alpha)\r
 136                 SetPPScore(PPSCORE_SPN);\r
 137 \r
 138         const unsigned uSeqCount1 = msa1.GetSeqCount();\r
 139         const unsigned uSeqCount2 = msa2.GetSeqCount();\r
 140         const unsigned uSumSeqCount = uSeqCount1 + uSeqCount2;\r
 141         MSA::SetIdCount(uSumSeqCount);\r
 142 \r
 143         ProfileProfile(msa1, msa2, msaOut);\r
 144 \r
 145         Progress("Writing output");\r
 146         MuscleOutput(msaOut);\r
 147         }\r