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