3 #include "objscore.h"
\r
14 static GAPINFO **g_Gaps;
\r
15 static GAPINFO *g_FreeList;
\r
16 static unsigned g_MaxSeqCount;
\r
17 static unsigned g_MaxColCount;
\r
18 static unsigned g_ColCount;
\r
19 static bool *g_ColDiff;
\r
21 static GAPINFO *NewGapInfo()
\r
23 if (0 == g_FreeList)
\r
25 const int NEWCOUNT = 256;
\r
26 GAPINFO *NewList = new GAPINFO[NEWCOUNT];
\r
27 g_FreeList = &NewList[0];
\r
28 for (int i = 0; i < NEWCOUNT-1; ++i)
\r
29 NewList[i].Next = &NewList[i+1];
\r
30 NewList[NEWCOUNT-1].Next = 0;
\r
32 GAPINFO *GI = g_FreeList;
\r
33 g_FreeList = g_FreeList->Next;
\r
37 static void FreeGapInfo(GAPINFO *GI)
\r
39 GI->Next = g_FreeList;
\r
43 // TODO: This could be much faster, no need to look
\r
45 static void FindIntersectingGaps(const MSA &msa, unsigned SeqIndex)
\r
47 const unsigned ColCount = msa.GetColCount();
\r
49 bool Intersects = false;
\r
50 unsigned Start = uInsane;
\r
51 for (unsigned Col = 0; Col <= ColCount; ++Col)
\r
53 bool Gap = ((Col != ColCount) && msa.IsGap(SeqIndex, Col));
\r
69 GAPINFO *GI = NewGapInfo();
\r
72 GI->Next = g_Gaps[SeqIndex];
\r
73 g_Gaps[SeqIndex] = GI;
\r
80 static SCORE Penalty(unsigned Length, bool Term)
\r
84 SCORE s1 = g_scoreGapOpen + g_scoreGapExtend*(Length - 1);
\r
86 SCORE s2 = g_scoreGapOpen2 + g_scoreGapExtend2*(Length - 1);
\r
95 //static SCORE ScorePair(unsigned Seq1, unsigned Seq2)
\r
99 // Log("ScorePair(%d,%d)\n", Seq1, Seq2);
\r
100 // Log("Gaps seq 1: ");
\r
101 // for (GAPINFO *GI = g_Gaps[Seq1]; GI; GI = GI->Next)
\r
102 // Log(" %d-%d", GI->Start, GI->End);
\r
104 // Log("Gaps seq 2: ");
\r
105 // for (GAPINFO *GI = g_Gaps[Seq2]; GI; GI = GI->Next)
\r
106 // Log(" %d-%d", GI->Start, GI->End);
\r
113 SCORE ScoreGaps(const MSA &msa, const unsigned DiffCols[], unsigned DiffColCount)
\r
117 Log("ScoreGaps\n");
\r
119 for (unsigned i = 0; i < DiffColCount; ++i)
\r
120 Log(" %u", DiffCols[i]);
\r
127 const unsigned SeqCount = msa.GetSeqCount();
\r
128 const unsigned ColCount = msa.GetColCount();
\r
129 g_ColCount = ColCount;
\r
131 if (SeqCount > g_MaxSeqCount)
\r
134 g_MaxSeqCount = SeqCount + 256;
\r
135 g_Gaps = new GAPINFO *[g_MaxSeqCount];
\r
137 memset(g_Gaps, 0, SeqCount*sizeof(GAPINFO *));
\r
139 if (ColCount > g_MaxColCount)
\r
141 delete[] g_ColDiff;
\r
142 g_MaxColCount = ColCount + 256;
\r
143 g_ColDiff = new bool[g_MaxColCount];
\r
146 memset(g_ColDiff, 0, g_ColCount*sizeof(bool));
\r
147 for (unsigned i = 0; i < DiffColCount; ++i)
\r
149 unsigned Col = DiffCols[i];
\r
150 assert(Col < ColCount);
\r
151 g_ColDiff[Col] = true;
\r
154 for (unsigned SeqIndex = 0; SeqIndex < SeqCount; ++SeqIndex)
\r
155 FindIntersectingGaps(msa, SeqIndex);
\r
160 Log("Intersecting gaps:\n");
\r
162 for (unsigned Col = 0; Col < ColCount; ++Col)
\r
163 Log("%c", g_ColDiff[Col] ? '*' : ' ');
\r
166 for (unsigned Col = 0; Col < ColCount; ++Col)
\r
169 for (unsigned Seq = 0; Seq < SeqCount; ++Seq)
\r
172 for (unsigned Col = 0; Col < ColCount; ++Col)
\r
173 Log("%c", msa.GetChar(Seq, Col));
\r
175 for (GAPINFO *GI = g_Gaps[Seq]; GI; GI = GI->Next)
\r
176 Log(" (%d,%d)", GI->Start, GI->End);
\r
177 Log(" >%s\n", msa.GetSeqName(Seq));
\r
184 for (unsigned Seq1 = 0; Seq1 < SeqCount; ++Seq1)
\r
186 const WEIGHT w1 = msa.GetSeqWeight(Seq1);
\r
187 for (unsigned Seq2 = Seq1 + 1; Seq2 < SeqCount; ++Seq2)
\r
189 const WEIGHT w2 = msa.GetSeqWeight(Seq2);
\r
190 // const SCORE Pair = ScorePair(Seq1, Seq2);
\r
191 const SCORE Pair = ScoreSeqPairGaps(msa, Seq1, msa, Seq2);
\r
192 Score += w1*w2*Pair;
\r
194 Log("Seq1=%u Seq2=%u ScorePair=%.4g w1=%.4g w2=%.4g Sum=%.4g\n",
\r
195 Seq1, Seq2, Pair, w1, w2, Score);
\r