3 #include "diaglist.h"
\r
7 const unsigned KTUP = 5;
\r
8 const unsigned KTUPS = 6*6*6*6*6;
\r
9 static unsigned TuplePos[KTUPS];
\r
11 static char *TupleToStr(int t)
\r
14 int t1, t2, t3, t4, t5;
\r
20 t5 = (t/(6*6*6*6))%6;
\r
30 static unsigned GetTuple(const ProfPos *PP, unsigned uPos)
\r
32 const unsigned t0 = PP[uPos].m_uResidueGroup;
\r
33 if (RESIDUE_GROUP_MULTIPLE == t0)
\r
36 const unsigned t1 = PP[uPos+1].m_uResidueGroup;
\r
37 if (RESIDUE_GROUP_MULTIPLE == t1)
\r
40 const unsigned t2 = PP[uPos+2].m_uResidueGroup;
\r
41 if (RESIDUE_GROUP_MULTIPLE == t2)
\r
44 const unsigned t3 = PP[uPos+3].m_uResidueGroup;
\r
45 if (RESIDUE_GROUP_MULTIPLE == t3)
\r
48 const unsigned t4 = PP[uPos+4].m_uResidueGroup;
\r
49 if (RESIDUE_GROUP_MULTIPLE == t4)
\r
52 return t0 + t1*6 + t2*6*6 + t3*6*6*6 + t4*6*6*6*6;
\r
55 void FindDiags(const ProfPos *PX, unsigned uLengthX, const ProfPos *PY,
\r
56 unsigned uLengthY, DiagList &DL)
\r
58 if (ALPHA_Amino != g_Alpha)
\r
59 Quit("FindDiags: requires amino acid alphabet");
\r
63 if (uLengthX < 12 || uLengthY < 12)
\r
66 // Set A to shorter profile, B to longer
\r
72 if (uLengthX < uLengthY)
\r
77 uLengthA = uLengthX;
\r
78 uLengthB = uLengthY;
\r
85 uLengthA = uLengthY;
\r
86 uLengthB = uLengthX;
\r
89 // Build tuple map for the longer profile, B
\r
90 if (uLengthB < KTUP)
\r
91 Quit("FindDiags: profile too short");
\r
93 memset(TuplePos, EMPTY, sizeof(TuplePos));
\r
95 for (unsigned uPos = 0; uPos < uLengthB - KTUP; ++uPos)
\r
97 const unsigned uTuple = GetTuple(PB, uPos);
\r
98 if (EMPTY == uTuple)
\r
100 TuplePos[uTuple] = uPos;
\r
104 for (unsigned uPosA = 0; uPosA < uLengthA - KTUP; ++uPosA)
\r
106 const unsigned uTuple = GetTuple(PA, uPosA);
\r
107 if (EMPTY == uTuple)
\r
109 const unsigned uPosB = TuplePos[uTuple];
\r
110 if (EMPTY == uPosB)
\r
113 // This tuple is found in both profiles
\r
114 unsigned uStartPosA = uPosA;
\r
115 unsigned uStartPosB = uPosB;
\r
117 // Try to extend the match forwards
\r
118 unsigned uEndPosA = uPosA + KTUP - 1;
\r
119 unsigned uEndPosB = uPosB + KTUP - 1;
\r
122 if (uLengthA - 1 == uEndPosA || uLengthB - 1 == uEndPosB)
\r
124 const unsigned uAAGroupA = PA[uEndPosA+1].m_uResidueGroup;
\r
125 if (RESIDUE_GROUP_MULTIPLE == uAAGroupA)
\r
127 const unsigned uAAGroupB = PB[uEndPosB+1].m_uResidueGroup;
\r
128 if (RESIDUE_GROUP_MULTIPLE == uAAGroupB)
\r
130 if (uAAGroupA != uAAGroupB)
\r
139 Log("Match: A %4u-%4u ", uStartPosA, uEndPosA);
\r
140 for (unsigned n = uStartPosA; n <= uEndPosA; ++n)
\r
141 Log("%c", 'A' + PA[n].m_uResidueGroup);
\r
143 Log(" B %4u-%4u ", uStartPosB, uEndPosB);
\r
144 for (unsigned n = uStartPosB; n <= uEndPosB; ++n)
\r
145 Log("%c", 'A' + PB[n].m_uResidueGroup);
\r
150 const unsigned uLength = uEndPosA - uStartPosA + 1;
\r
151 assert(uEndPosB - uStartPosB + 1 == uLength);
\r
153 if (uLength >= g_uMinDiagLength)
\r
156 DL.Add(uStartPosB, uStartPosA, uLength);
\r
158 DL.Add(uStartPosA, uStartPosB, uLength);
\r