7 // Textbook Smith-Waterman affine gap implementation.
\r
11 static const char *LocalScoreToStr(SCORE s)
\r
13 static char str[16];
\r
14 if (MINUS_INFINITY == s)
\r
16 sprintf(str, "%6.2f", s);
\r
20 static void ListDP(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,
\r
21 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
24 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
27 if (uPrefixLengthB > 0)
\r
28 c = ConsensusChar(PB[uPrefixLengthB - 1]);
\r
29 Log(" %4u:%c", uPrefixLengthB, c);
\r
32 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
35 if (uPrefixLengthA > 0)
\r
36 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
37 Log("%4u:%c ", uPrefixLengthA, c);
\r
38 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
39 Log(" %s", LocalScoreToStr(DPM(uPrefixLengthA, uPrefixLengthB)));
\r
44 SCORE SW(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
\r
45 unsigned uLengthB, PWPath &Path)
\r
47 assert(uLengthB > 0 && uLengthA > 0);
\r
49 const unsigned uPrefixCountA = uLengthA + 1;
\r
50 const unsigned uPrefixCountB = uLengthB + 1;
\r
52 // Allocate DP matrices
\r
53 const size_t LM = uPrefixCountA*uPrefixCountB;
\r
54 SCORE *DPM_ = new SCORE[LM];
\r
55 SCORE *DPD_ = new SCORE[LM];
\r
56 SCORE *DPI_ = new SCORE[LM];
\r
59 DPD(0, 0) = MINUS_INFINITY;
\r
60 DPI(0, 0) = MINUS_INFINITY;
\r
62 DPM(1, 0) = MINUS_INFINITY;
\r
63 DPD(1, 0) = MINUS_INFINITY;
\r
64 DPI(1, 0) = MINUS_INFINITY;
\r
66 DPM(0, 1) = MINUS_INFINITY;
\r
67 DPD(0, 1) = MINUS_INFINITY;
\r
68 DPI(0, 1) = MINUS_INFINITY;
\r
70 // Empty prefix of B is special case
\r
71 for (unsigned uPrefixLengthA = 2; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
73 // M=LetterA+LetterB, impossible with empty prefix
\r
74 DPM(uPrefixLengthA, 0) = MINUS_INFINITY;
\r
76 // D=LetterA+GapB, never optimal in local alignment with gap penalties
\r
77 DPD(uPrefixLengthA, 0) = MINUS_INFINITY;
\r
79 // I=GapA+LetterB, impossible with empty prefix
\r
80 DPI(uPrefixLengthA, 0) = MINUS_INFINITY;
\r
83 // Empty prefix of A is special case
\r
84 for (unsigned uPrefixLengthB = 2; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
86 // M=LetterA+LetterB, impossible with empty prefix
\r
87 DPM(0, uPrefixLengthB) = MINUS_INFINITY;
\r
89 // D=LetterA+GapB, impossible with empty prefix
\r
90 DPD(0, uPrefixLengthB) = MINUS_INFINITY;
\r
92 // I=GapA+LetterB, never optimal in local alignment with gap penalties
\r
93 DPI(0, uPrefixLengthB) = MINUS_INFINITY;
\r
96 SCORE scoreMax = MINUS_INFINITY;
\r
97 unsigned uPrefixLengthAMax = uInsane;
\r
98 unsigned uPrefixLengthBMax = uInsane;
\r
103 SCORE scoreGapCloseB = MINUS_INFINITY;
\r
104 for (unsigned uPrefixLengthB = 1; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
106 const ProfPos &PPB = PB[uPrefixLengthB - 1];
\r
108 SCORE scoreGapCloseA = MINUS_INFINITY;
\r
109 for (unsigned uPrefixLengthA = 1; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
111 const ProfPos &PPA = PA[uPrefixLengthA - 1];
\r
114 // Match M=LetterA+LetterB
\r
115 SCORE scoreLL = ScoreProfPos2(PPA, PPB);
\r
117 SCORE scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1);
\r
118 SCORE scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseA;
\r
119 SCORE scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseB;
\r
122 if (scoreMM >= scoreDM && scoreMM >= scoreIM)
\r
123 scoreBest = scoreMM;
\r
124 else if (scoreDM >= scoreMM && scoreDM >= scoreIM)
\r
125 scoreBest = scoreDM;
\r
128 assert(scoreIM >= scoreMM && scoreIM >= scoreDM);
\r
129 scoreBest = scoreIM;
\r
133 scoreBest += scoreLL;
\r
134 if (scoreBest > scoreMax)
\r
136 scoreMax = scoreBest;
\r
137 uPrefixLengthAMax = uPrefixLengthA;
\r
138 uPrefixLengthBMax = uPrefixLengthB;
\r
140 DPM(uPrefixLengthA, uPrefixLengthB) = scoreBest;
\r
144 // Delete D=LetterA+GapB
\r
145 SCORE scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) +
\r
146 PA[uPrefixLengthA-1].m_scoreGapOpen;
\r
147 SCORE scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB);
\r
150 if (scoreMD >= scoreDD)
\r
151 scoreBest = scoreMD;
\r
154 assert(scoreDD >= scoreMD);
\r
155 scoreBest = scoreDD;
\r
157 DPD(uPrefixLengthA, uPrefixLengthB) = scoreBest;
\r
160 // Insert I=GapA+LetterB
\r
162 SCORE scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) +
\r
163 PB[uPrefixLengthB - 1].m_scoreGapOpen;
\r
164 SCORE scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1);
\r
167 if (scoreMI >= scoreII)
\r
168 scoreBest = scoreMI;
\r
171 assert(scoreII > scoreMI);
\r
172 scoreBest = scoreII;
\r
174 DPI(uPrefixLengthA, uPrefixLengthB) = scoreBest;
\r
177 scoreGapCloseA = PPA.m_scoreGapClose;
\r
179 scoreGapCloseB = PPB.m_scoreGapClose;
\r
184 ListDP(DPM_, PA, PB, uPrefixLengthA, uPrefixLengthB);
\r
186 ListDP(DPD_, PA, PB, uPrefixLengthA, uPrefixLengthB);
\r
188 ListDP(DPI_, PA, PB, uPrefixLengthA, uPrefixLengthB);
\r
191 assert(scoreMax == DPM(uPrefixLengthAMax, uPrefixLengthBMax));
\r
192 TraceBackSW(PA, uLengthA, PB, uLengthB, DPM_, DPD_, DPI_,
\r
193 uPrefixLengthAMax, uPrefixLengthBMax, Path);
\r
196 SCORE scorePath = FastScorePath2(PA, uLengthA, PB, uLengthB, Path);
\r
198 Log("Score = %s Path = %s\n", LocalScoreToStr(scoreMax), LocalScoreToStr(scorePath));
\r