website/archive/binaries/mac/src/muscle/ppscore.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 g_bTracePPScore = false;\r
   9 MSA *g_ptrPPScoreMSA1 = 0;\r
  10 MSA *g_ptrPPScoreMSA2 = 0;\r
  11 \r
  12 static ProfPos *ProfileFromMSALocal(MSA &msa, Tree &tree)\r
  13         {\r
  14         const unsigned uSeqCount = msa.GetSeqCount();\r
  15         for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
  16                 msa.SetSeqId(uSeqIndex, uSeqIndex);\r
  17 \r
  18         TreeFromMSA(msa, tree, g_Cluster2, g_Distance2, g_Root1);\r
  19         SetMuscleTree(tree);\r
  20         return ProfileFromMSA(msa);\r
  21         }\r
  22 \r
  23 void PPScore()\r
  24         {\r
  25         if (0 == g_pstrFileName1 || 0 == g_pstrFileName2)\r
  26                 Quit("-ppscore needs -in1 and -in2");\r
  27 \r
  28         SetSeqWeightMethod(g_SeqWeight1);\r
  29 \r
  30         TextFile file1(g_pstrFileName1);\r
  31         TextFile file2(g_pstrFileName2);\r
  32 \r
  33         MSA msa1;\r
  34         MSA msa2;\r
  35 \r
  36         msa1.FromFile(file1);\r
  37         msa2.FromFile(file2);\r
  38 \r
  39         const unsigned uLength1 = msa1.GetColCount();\r
  40         const unsigned uLength2 = msa2.GetColCount();\r
  41 \r
  42         if (uLength1 != uLength2)\r
  43                 Quit("Profiles must have the same length");\r
  44 \r
  45         ALPHA Alpha = ALPHA_Undefined;\r
  46         switch (g_SeqType)\r
  47                 {\r
  48         case SEQTYPE_Auto:\r
  49                 Alpha = msa1.GuessAlpha();\r
  50                 break;\r
  51 \r
  52         case SEQTYPE_Protein:\r
  53                 Alpha = ALPHA_Amino;\r
  54                 break;\r
  55 \r
  56         case SEQTYPE_DNA:\r
  57                 Alpha = ALPHA_DNA;\r
  58                 break;\r
  59 \r
  60         case SEQTYPE_RNA:\r
  61                 Alpha = ALPHA_RNA;\r
  62                 break;\r
  63 \r
  64         default:\r
  65                 Quit("Invalid SeqType");\r
  66                 }\r
  67         SetAlpha(Alpha);\r
  68 \r
  69         msa1.FixAlpha();\r
  70         msa2.FixAlpha();\r
  71 \r
  72         if (ALPHA_DNA == Alpha || ALPHA_RNA == Alpha)\r
  73                 SetPPScore(PPSCORE_SPN);\r
  74 \r
  75         const unsigned uSeqCount1 = msa1.GetSeqCount();\r
  76         const unsigned uSeqCount2 = msa2.GetSeqCount();\r
  77         const unsigned uMaxSeqCount = (uSeqCount1 > uSeqCount2 ? uSeqCount1 : uSeqCount2);\r
  78         MSA::SetIdCount(uMaxSeqCount);\r
  79 \r
  80         Tree tree1;\r
  81         Tree tree2;\r
  82         ProfPos *Prof1 = ProfileFromMSALocal(msa1, tree1);\r
  83         ProfPos *Prof2 = ProfileFromMSALocal(msa2, tree2);\r
  84 \r
  85         g_bTracePPScore = true;\r
  86         g_ptrPPScoreMSA1 = &msa1;\r
  87         g_ptrPPScoreMSA2 = &msa2;\r
  88 \r
  89         SCORE Score = ObjScoreDP_Profs(Prof1, Prof2, uLength1);\r
  90 \r
  91         Log("Score=%.4g\n", Score);\r
  92         printf("Score=%.4g\n", Score);\r
  93         }\r