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
|
#include "libMUSCLE/muscle.h"
#include "libMUSCLE/tree.h"
#include "libMUSCLE/profile.h"
#include "libMUSCLE/msa.h"
#include "libMUSCLE/seqvect.h"
#include "libMUSCLE/pwpath.h"
namespace muscle {
static void DoSeq(Seq &s, unsigned uSeqIndex, const ProfPos *RootProf,
unsigned uRootProfLength, MSA &msaOut)
{
MSA msaSeq;
msaSeq.FromSeq(s);
const unsigned uSeqLength = s.Length();
MSA msaDummy;
msaDummy.SetSize(1, uRootProfLength);
msaDummy.SetSeqId(0, 0);
msaDummy.SetSeqName(0, "Dummy0");
for (unsigned uColIndex = 0; uColIndex < uRootProfLength; ++uColIndex)
msaDummy.SetChar(0, uColIndex, '?');
ProfPos *SeqProf = ProfileFromMSA(msaSeq);
for (unsigned uColIndex = 0; uColIndex < uSeqLength; ++uColIndex)
{
ProfPos &PP = SeqProf[uColIndex];
PP.m_scoreGapOpen = MINUS_INFINITY;
PP.m_scoreGapClose = MINUS_INFINITY;
}
ProfPos *ProfOut;
unsigned uLengthOut;
PWPath Path;
AlignTwoProfs(SeqProf, uSeqLength, 1.0, RootProf, uRootProfLength, 1.0,
Path, &ProfOut, &uLengthOut);
assert(uLengthOut = uRootProfLength);
delete[] ProfOut;
MSA msaCombined;
AlignTwoMSAsGivenPath(Path, msaSeq, msaDummy, msaCombined);
msaCombined.LogMe();
msaOut.SetSeqName(uSeqIndex, s.GetName());
msaOut.SetSeqId(uSeqIndex, s.GetId());
for (unsigned uColIndex = 0; uColIndex < uRootProfLength; ++uColIndex)
msaOut.SetChar(uSeqIndex, uColIndex, msaCombined.GetChar(0, uColIndex));
}
// Steven Brenner's O(NL^2) proposal for creating a root alignment
// Align each sequence to the profile at the root.
// Compare the e-string solution, which is O(NL log N).
void MakeRootMSABrenner(SeqVect &v, const Tree &GuideTree, ProgNode Nodes[],
MSA &a)
{
const unsigned uSeqCount = v.Length();
const unsigned uRootNodeIndex = GuideTree.GetRootNodeIndex();
const ProfPos *RootProfile = Nodes[uRootNodeIndex].m_Prof;
const unsigned uRootColCount = Nodes[uRootNodeIndex].m_uLength;
a.SetSize(uSeqCount, uRootColCount);
for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)
DoSeq(*v[uSeqIndex], uSeqIndex, RootProfile, uRootColCount, a);
}
}
|
