website/archive/binaries/mac/src/muscle/finddiags.cpp

   1 #include "muscle.h"\r
   2 #include "profile.h"\r
   3 #include "diaglist.h"\r
   4 \r
   5 #define TRACE   0\r
   6 \r
   7 const unsigned KTUP = 5;\r
   8 const unsigned KTUPS = 6*6*6*6*6;\r
   9 static unsigned TuplePos[KTUPS];\r
  10 \r
  11 static char *TupleToStr(int t)\r
  12         {\r
  13         static char s[7];\r
  14         int t1, t2, t3, t4, t5;\r
  15 \r
  16         t1 = t%6;\r
  17         t2 = (t/6)%6;\r
  18         t3 = (t/(6*6))%6;\r
  19         t4 = (t/(6*6*6))%6;\r
  20         t5 = (t/(6*6*6*6))%6;\r
  21 \r
  22         s[4] = '0' + t1;\r
  23         s[3] = '0' + t2;\r
  24         s[2] = '0' + t3;\r
  25         s[1] = '0' + t4;\r
  26         s[0] = '0' + t5;\r
  27         return s;\r
  28         }\r
  29 \r
  30 static unsigned GetTuple(const ProfPos *PP, unsigned uPos)\r
  31         {\r
  32         const unsigned t0 = PP[uPos].m_uResidueGroup;\r
  33         if (RESIDUE_GROUP_MULTIPLE == t0)\r
  34                 return EMPTY;\r
  35 \r
  36         const unsigned t1 = PP[uPos+1].m_uResidueGroup;\r
  37         if (RESIDUE_GROUP_MULTIPLE == t1)\r
  38                 return EMPTY;\r
  39 \r
  40         const unsigned t2 = PP[uPos+2].m_uResidueGroup;\r
  41         if (RESIDUE_GROUP_MULTIPLE == t2)\r
  42                 return EMPTY;\r
  43 \r
  44         const unsigned t3 = PP[uPos+3].m_uResidueGroup;\r
  45         if (RESIDUE_GROUP_MULTIPLE == t3)\r
  46                 return EMPTY;\r
  47 \r
  48         const unsigned t4 = PP[uPos+4].m_uResidueGroup;\r
  49         if (RESIDUE_GROUP_MULTIPLE == t4)\r
  50                 return EMPTY;\r
  51 \r
  52         return t0 + t1*6 + t2*6*6 + t3*6*6*6 + t4*6*6*6*6;\r
  53         }\r
  54 \r
  55 void FindDiags(const ProfPos *PX, unsigned uLengthX, const ProfPos *PY,\r
  56   unsigned uLengthY, DiagList &DL)\r
  57         {\r
  58         if (ALPHA_Amino != g_Alpha)\r
  59                 Quit("FindDiags: requires amino acid alphabet");\r
  60 \r
  61         DL.Clear();\r
  62 \r
  63         if (uLengthX < 12 || uLengthY < 12)\r
  64                 return;\r
  65 \r
  66 // Set A to shorter profile, B to longer\r
  67         const ProfPos *PA;\r
  68         const ProfPos *PB;\r
  69         unsigned uLengthA;\r
  70         unsigned uLengthB;\r
  71         bool bSwap;\r
  72         if (uLengthX < uLengthY)\r
  73                 {\r
  74                 bSwap = false;\r
  75                 PA = PX;\r
  76                 PB = PY;\r
  77                 uLengthA = uLengthX;\r
  78                 uLengthB = uLengthY;\r
  79                 }\r
  80         else\r
  81                 {\r
  82                 bSwap = true;\r
  83                 PA = PY;\r
  84                 PB = PX;\r
  85                 uLengthA = uLengthY;\r
  86                 uLengthB = uLengthX;\r
  87                 }\r
  88 \r
  89 // Build tuple map for the longer profile, B\r
  90         if (uLengthB < KTUP)\r
  91                 Quit("FindDiags: profile too short");\r
  92 \r
  93         memset(TuplePos, EMPTY, sizeof(TuplePos));\r
  94 \r
  95         for (unsigned uPos = 0; uPos < uLengthB - KTUP; ++uPos)\r
  96                 {\r
  97                 const unsigned uTuple = GetTuple(PB, uPos);\r
  98                 if (EMPTY == uTuple)\r
  99                         continue;\r
 100                 TuplePos[uTuple] = uPos;\r
 101                 }\r
 102 \r
 103 // Find matches\r
 104         for (unsigned uPosA = 0; uPosA < uLengthA - KTUP; ++uPosA)\r
 105                 {\r
 106                 const unsigned uTuple = GetTuple(PA, uPosA);\r
 107                 if (EMPTY == uTuple)\r
 108                         continue;\r
 109                 const unsigned uPosB = TuplePos[uTuple];\r
 110                 if (EMPTY == uPosB)\r
 111                         continue;\r
 112 \r
 113         // This tuple is found in both profiles\r
 114                 unsigned uStartPosA = uPosA;\r
 115                 unsigned uStartPosB = uPosB;\r
 116 \r
 117         // Try to extend the match forwards\r
 118                 unsigned uEndPosA = uPosA + KTUP - 1;\r
 119                 unsigned uEndPosB = uPosB + KTUP - 1;\r
 120                 for (;;)\r
 121                         {\r
 122                         if (uLengthA - 1 == uEndPosA || uLengthB - 1 == uEndPosB)\r
 123                                 break;\r
 124                         const unsigned uAAGroupA = PA[uEndPosA+1].m_uResidueGroup;\r
 125                         if (RESIDUE_GROUP_MULTIPLE == uAAGroupA)\r
 126                                 break;\r
 127                         const unsigned uAAGroupB = PB[uEndPosB+1].m_uResidueGroup;\r
 128                         if (RESIDUE_GROUP_MULTIPLE == uAAGroupB)\r
 129                                 break;\r
 130                         if (uAAGroupA != uAAGroupB)\r
 131                                 break;\r
 132                         ++uEndPosA;\r
 133                         ++uEndPosB;\r
 134                         }\r
 135                 uPosA = uEndPosA;\r
 136 \r
 137 #if     TRACE\r
 138                 {\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
 142                 Log("\n");\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
 146                 Log("\n");\r
 147                 }\r
 148 #endif\r
 149 \r
 150                 const unsigned uLength = uEndPosA - uStartPosA + 1;\r
 151                 assert(uEndPosB - uStartPosB + 1 == uLength);\r
 152 \r
 153                 if (uLength >= g_uMinDiagLength)\r
 154                         {\r
 155                         if (bSwap)\r
 156                                 DL.Add(uStartPosB, uStartPosA, uLength);\r
 157                         else\r
 158                                 DL.Add(uStartPosA, uStartPosB, uLength);\r
 159                         }\r
 160                 }\r
 161         }\r