1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
|
#include "muscle.h"
#include "msa.h"
#include "profile.h"
#include "pwpath.h"
#include "textfile.h"
#include "timing.h"
SCORE GlobalAlign4(ProfPos *PA, unsigned uLengthA, ProfPos *PB,
unsigned uLengthB, PWPath &Path);
static SCORE GlobalAlignXX(int NPenalty, int CPenalty, ProfPos *PA, unsigned uLengthA,
ProfPos *PB, unsigned uLengthB, PWPath &Path, SCORE *ptrMatch);
static SCORE MatchScore(ProfPos *PA, unsigned uLengthA, ProfPos *PB,
unsigned uLengthB, PWPath &Path);
SCORE AlignTwoMSAs(const MSA &msa1, const MSA &msa2, MSA &msaOut, PWPath &Path,
bool bLockLeft, bool bLockRight)
{
const unsigned uLengthA = msa1.GetColCount();
const unsigned uLengthB = msa2.GetColCount();
SCORE scoreGapOpen = (SCORE) g_scoreGapOpen;
bool bTermGapsFreeA = TermGapsFree(uLengthA, uLengthB);
bool bTermGapsFreeB = TermGapsFree(uLengthB, uLengthA);
ProfPos *PA = ProfileFromMSA(msa1, scoreGapOpen, bTermGapsFreeA);
ProfPos *PB = ProfileFromMSA(msa2, scoreGapOpen, bTermGapsFreeB);
if (bLockLeft)
{
PA[0].m_scoreGapOpen = MINUS_INFINITY;
PB[0].m_scoreGapOpen = MINUS_INFINITY;
}
if (bLockRight)
{
PA[uLengthA-1].m_scoreGapClose = MINUS_INFINITY;
PB[uLengthB-1].m_scoreGapClose = MINUS_INFINITY;
}
float r = (float) uLengthA/ (float) (uLengthB + 1); // +1 to prevent div 0
if (r < 1)
r = 1/r;
SCORE Score;
if (g_bTermGaps4 && !bLockLeft && !bLockRight && r > 1.2)
Score = GlobalAlign4(PA, uLengthA, PB, uLengthB, Path);
else
Score = GlobalAlign(PA, uLengthA, PB, uLengthB, Path);
AlignTwoMSAsGivenPath(Path, msa1, msa2, msaOut);
delete[] PA;
delete[] PB;
return Score;
}
SCORE GlobalAlign4(ProfPos *PA, unsigned uLengthA, ProfPos *PB,
unsigned uLengthB, PWPath &Path)
{
SCORE Score11;
SCORE Score10;
SCORE Score01;
SCORE Score00;
SCORE Match11;
SCORE Match10;
SCORE Match01;
SCORE Match00;
PWPath Path11;
PWPath Path10;
PWPath Path01;
PWPath Path00;
Score11 = GlobalAlignXX(1, 1, PA, uLengthA, PB, uLengthB, Path11, &Match11);
Score10 = GlobalAlignXX(1, 0, PA, uLengthA, PB, uLengthB, Path10, &Match10);
Score01 = GlobalAlignXX(0, 1, PA, uLengthA, PB, uLengthB, Path01, &Match01);
Score00 = GlobalAlignXX(0, 0, PA, uLengthA, PB, uLengthB, Path00, &Match00);
SCORE MaxMatch = Match11;
int Best = 1;
if (Match10 > MaxMatch)
{
Best = 2;
MaxMatch = Match10;
}
if (Match01 > MaxMatch)
{
Best = 3;
MaxMatch = Match01;
}
if (Match00 > MaxMatch)
{
Best = 4;
MaxMatch = Match00;
}
switch (Best)
{
case 1:
Path.Copy(Path11);
return Score11;
case 2:
Path.Copy(Path10);
return Score10;
case 3:
Path.Copy(Path01);
return Score01;
case 4:
Path.Copy(Path00);
return Score00;
}
Quit("GlobalAlign4 invalid Best");
return MINUS_INFINITY;
}
static SCORE GlobalAlignXX(int NPenalty, int CPenalty, ProfPos *PA,
unsigned uLengthA, ProfPos *PB, unsigned uLengthB, PWPath &Path, SCORE *ptrMatch)
{
bool bALonger = (uLengthA > uLengthB);
ProfPos &PPN = (bALonger ? PA[0] : PB[0]);
ProfPos &PPC = (bALonger ? PA[uLengthA - 1] : PB[uLengthB - 1]);
if (NPenalty)
PPN.m_scoreGapOpen = g_scoreGapOpen;
else
PPN.m_scoreGapOpen = 0;
if (CPenalty)
PPC.m_scoreGapClose = g_scoreGapOpen;
else
PPC.m_scoreGapClose = 0;
SCORE Score = GlobalAlign(PA, uLengthA, PB, uLengthB, Path);
*ptrMatch = MatchScore(PA, uLengthA, PB, uLengthB, Path);
return Score;
}
static SCORE MatchScore(ProfPos *PA, unsigned uLengthA, ProfPos *PB,
unsigned uLengthB, PWPath &Path)
{
SCORE Score = 0;
const unsigned uEdgeCount = Path.GetEdgeCount();
for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)
{
const PWEdge &Edge = Path.GetEdge(uEdgeIndex);
if (Edge.cType != 'M')
continue;
const ProfPos &PPA = PA[Edge.uPrefixLengthA - 1];
const ProfPos &PPB = PB[Edge.uPrefixLengthB - 1];
Score += ScoreProfPos2(PPA, PPB);
}
return Score;
}
|
