4 #include "objscore.h"
\r
9 #define TEST_SPFAST 0
\r
11 static SCORE GapPenalty(unsigned uLength, bool Term, SCORE g, SCORE e)
\r
15 // switch (g_TermGap)
\r
17 // case TERMGAP_Full:
\r
18 // return g + (uLength - 1)*e;
\r
20 // case TERMGAP_Half:
\r
21 // return g/2 + (uLength - 1)*e;
\r
23 // case TERMGAP_Ext:
\r
24 // return uLength*e;
\r
26 // Quit("Bad termgap");
\r
29 // return g + (uLength - 1)*e;
\r
30 //return MINUS_INFINITY;
\r
31 return g + (uLength - 1)*e;
\r
34 static SCORE GapPenalty(unsigned uLength, bool Term)
\r
36 SCORE s1 = GapPenalty(uLength, Term, g_scoreGapOpen, g_scoreGapExtend);
\r
38 SCORE s2 = GapPenalty(uLength, Term, g_scoreGapOpen2, g_scoreGapExtend2);
\r
47 static const MSA *g_ptrMSA1;
\r
48 static const MSA *g_ptrMSA2;
\r
49 static unsigned g_uSeqIndex1;
\r
50 static unsigned g_uSeqIndex2;
\r
52 static void LogGap(unsigned uStart, unsigned uEnd, unsigned uGapLength,
\r
53 bool bNTerm, bool bCTerm)
\r
55 Log("%16.16s ", "");
\r
56 for (unsigned i = 0; i < uStart; ++i)
\r
58 unsigned uMyLength = 0;
\r
59 for (unsigned i = uStart; i <= uEnd; ++i)
\r
61 bool bGap1 = g_ptrMSA1->IsGap(g_uSeqIndex1, i);
\r
62 bool bGap2 = g_ptrMSA2->IsGap(g_uSeqIndex2, i);
\r
63 if (!bGap1 && !bGap2)
\r
64 Quit("Error -- neither gapping");
\r
73 SCORE s = GapPenalty(uGapLength, bNTerm || bCTerm);
\r
74 Log(" L=%d N%d C%d s=%.3g", uGapLength, bNTerm, bCTerm, s);
\r
76 if (uMyLength != uGapLength)
\r
77 Quit("Lengths differ");
\r
81 static SCORE ScoreSeqPair(const MSA &msa1, unsigned uSeqIndex1,
\r
82 const MSA &msa2, unsigned uSeqIndex2, SCORE *ptrLetters, SCORE *ptrGaps)
\r
86 g_uSeqIndex1 = uSeqIndex1;
\r
87 g_uSeqIndex2 = uSeqIndex2;
\r
89 const unsigned uColCount = msa1.GetColCount();
\r
90 const unsigned uColCount2 = msa2.GetColCount();
\r
91 if (uColCount != uColCount2)
\r
92 Quit("ScoreSeqPair, different lengths");
\r
95 Log("ScoreSeqPair\n");
\r
96 Log("%16.16s ", msa1.GetSeqName(uSeqIndex1));
\r
97 for (unsigned i = 0; i < uColCount; ++i)
\r
98 Log("%c", msa1.GetChar(uSeqIndex1, i));
\r
100 Log("%16.16s ", msa2.GetSeqName(uSeqIndex2));
\r
101 for (unsigned i = 0; i < uColCount; ++i)
\r
102 Log("%c", msa1.GetChar(uSeqIndex2, i));
\r
106 SCORE scoreTotal = 0;
\r
108 // Substitution scores
\r
109 unsigned uFirstLetter1 = uInsane;
\r
110 unsigned uFirstLetter2 = uInsane;
\r
111 unsigned uLastLetter1 = uInsane;
\r
112 unsigned uLastLetter2 = uInsane;
\r
113 for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)
\r
115 bool bGap1 = msa1.IsGap(uSeqIndex1, uColIndex);
\r
116 bool bGap2 = msa2.IsGap(uSeqIndex2, uColIndex);
\r
117 bool bWildcard1 = msa1.IsWildcard(uSeqIndex1, uColIndex);
\r
118 bool bWildcard2 = msa2.IsWildcard(uSeqIndex2, uColIndex);
\r
122 if (uInsane == uFirstLetter1)
\r
123 uFirstLetter1 = uColIndex;
\r
124 uLastLetter1 = uColIndex;
\r
128 if (uInsane == uFirstLetter2)
\r
129 uFirstLetter2 = uColIndex;
\r
130 uLastLetter2 = uColIndex;
\r
133 if (bGap1 || bGap2 || bWildcard1 || bWildcard2)
\r
136 unsigned uLetter1 = msa1.GetLetter(uSeqIndex1, uColIndex);
\r
137 unsigned uLetter2 = msa2.GetLetter(uSeqIndex2, uColIndex);
\r
139 SCORE scoreMatch = (*g_ptrScoreMatrix)[uLetter1][uLetter2];
\r
140 scoreTotal += scoreMatch;
\r
142 Log("%c <-> %c = %7.1f %10.1f\n",
\r
143 msa1.GetChar(uSeqIndex1, uColIndex),
\r
144 msa2.GetChar(uSeqIndex2, uColIndex),
\r
150 *ptrLetters = scoreTotal;
\r
153 unsigned uGapLength = uInsane;
\r
154 unsigned uGapStartCol = uInsane;
\r
155 bool bGapping1 = false;
\r
156 bool bGapping2 = false;
\r
158 for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)
\r
160 bool bGap1 = msa1.IsGap(uSeqIndex1, uColIndex);
\r
161 bool bGap2 = msa2.IsGap(uSeqIndex2, uColIndex);
\r
163 if (bGap1 && bGap2)
\r
173 bool bNTerm = (uFirstLetter2 == uGapStartCol);
\r
174 bool bCTerm = (uLastLetter2 + 1 == uColIndex);
\r
175 SCORE scoreGap = GapPenalty(uGapLength, bNTerm || bCTerm);
\r
176 scoreTotal += scoreGap;
\r
178 LogGap(uGapStartCol, uColIndex - 1, uGapLength, bNTerm, bCTerm);
\r
179 Log("GAP %7.1f %10.1f\n",
\r
190 uGapStartCol = uColIndex;
\r
204 bool bNTerm = (uFirstLetter1 == uGapStartCol);
\r
205 bool bCTerm = (uLastLetter1 + 1 == uColIndex);
\r
206 SCORE scoreGap = GapPenalty(uGapLength, bNTerm || bCTerm);
\r
207 scoreTotal += scoreGap;
\r
209 LogGap(uGapStartCol, uColIndex - 1, uGapLength, bNTerm, bCTerm);
\r
210 Log("GAP %7.1f %10.1f\n",
\r
220 uGapStartCol = uColIndex;
\r
227 if (bGapping1 || bGapping2)
\r
229 SCORE scoreGap = GapPenalty(uGapLength, true);
\r
230 scoreTotal += scoreGap;
\r
232 LogGap(uGapStartCol, uColCount - 1, uGapLength, false, true);
\r
233 Log("GAP %7.1f %10.1f\n",
\r
238 *ptrGaps = scoreTotal - *ptrLetters;
\r
242 // The usual sum-of-pairs objective score: sum the score
\r
243 // of the alignment of each pair of sequences.
\r
244 SCORE ObjScoreDA(const MSA &msa, SCORE *ptrLetters, SCORE *ptrGaps)
\r
246 const unsigned uSeqCount = msa.GetSeqCount();
\r
247 SCORE scoreTotal = 0;
\r
248 unsigned uPairCount = 0;
\r
251 Log(" Score Weight Weight Total\n");
\r
252 Log("---------- ------ ------ ----------\n");
\r
254 SCORE TotalLetters = 0;
\r
255 SCORE TotalGaps = 0;
\r
256 for (unsigned uSeqIndex1 = 0; uSeqIndex1 < uSeqCount; ++uSeqIndex1)
\r
258 const WEIGHT w1 = msa.GetSeqWeight(uSeqIndex1);
\r
259 for (unsigned uSeqIndex2 = uSeqIndex1 + 1; uSeqIndex2 < uSeqCount; ++uSeqIndex2)
\r
261 const WEIGHT w2 = msa.GetSeqWeight(uSeqIndex2);
\r
262 const WEIGHT w = w1*w2;
\r
265 SCORE scorePair = ScoreSeqPair(msa, uSeqIndex1, msa, uSeqIndex2,
\r
267 scoreTotal += w1*w2*scorePair;
\r
268 TotalLetters += w1*w2*Letters;
\r
269 TotalGaps += w1*w2*Gaps;
\r
272 Log("%10.2f %6.3f %6.3f %10.2f %d=%s %d=%s\n",
\r
278 msa.GetSeqName(uSeqIndex1),
\r
280 msa.GetSeqName(uSeqIndex2));
\r
284 *ptrLetters = TotalLetters;
\r
285 *ptrGaps = TotalGaps;
\r
289 #endif // DOUBLE_AFFINE
\r