binaries/src/muscle/sptest.cpp

   1 #include "muscle.h"\r
   2 #include "objscore.h"\r
   3 #include "msa.h"\r
   4 #include "textfile.h"\r
   5 #include "pwpath.h"\r
   6 \r
   7 const unsigned INDELS = 1;\r
   8 \r
   9 static void GetPos(const char Str[], unsigned L, int *pi1, int *pi2)\r
  10         {\r
  11         int i1;\r
  12         for (;;)\r
  13                 {\r
  14                 i1 = rand()%(L-2) + 1;\r
  15                 if (Str[i1] == 'M')\r
  16                         break;\r
  17                 }\r
  18         int i2;\r
  19         for (;;)\r
  20                 {\r
  21                 i2 = rand()%(L-2) + 1;\r
  22                 if (i1 != i2 && Str[i2] == 'M')\r
  23                         break;\r
  24                 }\r
  25         *pi1 = i1;\r
  26         *pi2 = i2;\r
  27         }\r
  28 \r
  29 static void MakePath(unsigned uSeqLength, unsigned uIndelCount, char Str[])\r
  30         {\r
  31         unsigned uPathLength = uSeqLength + uIndelCount;\r
  32         for (unsigned i = 0; i < uPathLength; ++i)\r
  33                 Str[i] = 'M';\r
  34 \r
  35         for (unsigned i = 0; i < uIndelCount; ++i)\r
  36                 {\r
  37                 int i1, i2;\r
  38                 GetPos(Str, uPathLength, &i1, &i2);\r
  39                 Str[i1] = 'D';\r
  40                 Str[i2] = 'I';\r
  41                 }\r
  42 \r
  43         Str[uPathLength] = 0;\r
  44         Log("MakePath=%s\n", Str);\r
  45         }\r
  46 \r
  47 void SPTest()\r
  48         {\r
  49         SetPPScore(PPSCORE_SV);\r
  50 \r
  51         SetListFileName("c:\\tmp\\muscle.log", false);\r
  52 \r
  53         TextFile file1("c:\\tmp\\msa1.afa");\r
  54         TextFile file2("c:\\tmp\\msa2.afa");\r
  55 \r
  56         MSA msa1;\r
  57         MSA msa2;\r
  58 \r
  59         msa1.FromFile(file1);\r
  60         msa2.FromFile(file2);\r
  61 \r
  62         Log("msa1=\n");\r
  63         msa1.LogMe();\r
  64         Log("msa2=\n");\r
  65         msa2.LogMe();\r
  66 \r
  67         const unsigned uColCount = msa1.GetColCount();\r
  68         if (msa2.GetColCount() != uColCount)\r
  69                 Quit("Different lengths");\r
  70 \r
  71         const unsigned uSeqCount1 = msa1.GetSeqCount();\r
  72         const unsigned uSeqCount2 = msa2.GetSeqCount();\r
  73         const unsigned uSeqCount = uSeqCount1 + uSeqCount2;\r
  74 \r
  75         MSA::SetIdCount(uSeqCount);\r
  76 \r
  77         for (unsigned uSeqIndex1 = 0; uSeqIndex1 < uSeqCount1; ++uSeqIndex1)\r
  78                 {\r
  79                 msa1.SetSeqWeight(uSeqIndex1, 1.0);\r
  80                 msa1.SetSeqId(uSeqIndex1, uSeqIndex1);\r
  81                 }\r
  82 \r
  83         for (unsigned uSeqIndex2 = 0; uSeqIndex2 < uSeqCount2; ++uSeqIndex2)\r
  84                 {\r
  85                 msa2.SetSeqWeight(uSeqIndex2, 1.0);\r
  86                 msa2.SetSeqId(uSeqIndex2, uSeqCount1 + uSeqIndex2);\r
  87                 }\r
  88 \r
  89         MSA alnA;\r
  90         MSA alnB;\r
  91 \r
  92         char strPathA[1024];\r
  93         char strPathB[1024];\r
  94         MakePath(uColCount, INDELS, strPathA);\r
  95         MakePath(uColCount, INDELS, strPathB);\r
  96 \r
  97         PWPath PathA;\r
  98         PWPath PathB;\r
  99         PathA.FromStr(strPathA);\r
 100         PathB.FromStr(strPathB);\r
 101 \r
 102         Log("PathA=\n");\r
 103         PathA.LogMe();\r
 104         Log("PathB=\n");\r
 105         PathB.LogMe();\r
 106 \r
 107         AlignTwoMSAsGivenPath(PathA, msa1, msa2, alnA);\r
 108         AlignTwoMSAsGivenPath(PathB, msa1, msa2, alnB);\r
 109 \r
 110         for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
 111                 {\r
 112                 alnA.SetSeqWeight(uSeqIndex, 1.0);\r
 113                 alnB.SetSeqWeight(uSeqIndex, 1.0);\r
 114                 }\r
 115 \r
 116         unsigned Seqs1[1024];\r
 117         unsigned Seqs2[1024];\r
 118 \r
 119         for (unsigned uSeqIndex1 = 0; uSeqIndex1 < uSeqCount1; ++uSeqIndex1)\r
 120                 Seqs1[uSeqIndex1] = uSeqIndex1;\r
 121 \r
 122         for (unsigned uSeqIndex2 = 0; uSeqIndex2 < uSeqCount2; ++uSeqIndex2)\r
 123                 Seqs2[uSeqIndex2] = uSeqCount1 + uSeqIndex2;\r
 124 \r
 125         MSA msaA1;\r
 126         MSA msaA2;\r
 127         MSA msaB1;\r
 128         MSA msaB2;\r
 129         MSAFromSeqSubset(alnA, Seqs1, uSeqCount1, msaA1);\r
 130         MSAFromSeqSubset(alnB, Seqs1, uSeqCount1, msaB1);\r
 131         MSAFromSeqSubset(alnA, Seqs2, uSeqCount2, msaA2);\r
 132         MSAFromSeqSubset(alnB, Seqs2, uSeqCount2, msaB2);\r
 133 \r
 134         for (unsigned uSeqIndex1 = 0; uSeqIndex1 < uSeqCount1; ++uSeqIndex1)\r
 135                 {\r
 136                 msaA1.SetSeqWeight(uSeqIndex1, 1.0);\r
 137                 msaB1.SetSeqWeight(uSeqIndex1, 1.0);\r
 138                 }\r
 139 \r
 140         for (unsigned uSeqIndex2 = 0; uSeqIndex2 < uSeqCount2; ++uSeqIndex2)\r
 141                 {\r
 142                 msaA2.SetSeqWeight(uSeqIndex2, 1.0);\r
 143                 msaB2.SetSeqWeight(uSeqIndex2, 1.0);\r
 144                 }\r
 145 \r
 146         Log("msaA1=\n");\r
 147         msaA1.LogMe();\r
 148 \r
 149         Log("msaB1=\n");\r
 150         msaB1.LogMe();\r
 151 \r
 152         Log("msaA2=\n");\r
 153         msaA2.LogMe();\r
 154 \r
 155         Log("msaB2=\n");\r
 156         msaB2.LogMe();\r
 157 \r
 158         Log("alnA=\n");\r
 159         alnA.LogMe();\r
 160 \r
 161         Log("AlnB=\n");\r
 162         alnB.LogMe();\r
 163 \r
 164         Log("\nSPA\n---\n");\r
 165         SCORE SPA = ObjScoreSP(alnA);\r
 166         Log("\nSPB\n---\n");\r
 167         SCORE SPB = ObjScoreSP(alnB);\r
 168 \r
 169         Log("\nXPA\n---\n");\r
 170         SCORE XPA = ObjScoreXP(msaA1, msaA2);\r
 171         Log("\nXPB\n---\n");\r
 172         SCORE XPB = ObjScoreXP(msaB1, msaB2);\r
 173 \r
 174         Log("SPA=%.4g SPB=%.4g Diff=%.4g\n", SPA, SPB, SPA - SPB);\r
 175         Log("XPA=%.4g XPB=%.4g Diff=%.4g\n", XPA, XPB, XPA - XPB);\r
 176         }\r