--- /dev/null
+#include "muscle.h"\r
+#include "distfunc.h"\r
+#include "seqvect.h"\r
+#include <math.h>\r
+\r
+const unsigned TRIPLE_COUNT = 20*20*20;\r
+\r
+struct TripleCount\r
+ {\r
+ unsigned m_uSeqCount; // How many sequences have this triple?\r
+ unsigned short *m_Counts; // m_Counts[s] = nr of times triple found in seq s\r
+ };\r
+static TripleCount *TripleCounts;\r
+\r
+// WARNING: Sequences MUST be stripped of gaps and upper case!\r
+void DistKmer20_3(const SeqVect &v, DistFunc &DF)\r
+ {\r
+ const unsigned uSeqCount = v.Length();\r
+\r
+ DF.SetCount(uSeqCount);\r
+ if (0 == uSeqCount)\r
+ return;\r
+ for (unsigned uSeq1 = 0; uSeq1 < uSeqCount; ++uSeq1)\r
+ {\r
+ DF.SetDist(uSeq1, uSeq1, 0);\r
+ for (unsigned uSeq2 = 0; uSeq2 < uSeq1; ++uSeq2)\r
+ DF.SetDist(uSeq1, uSeq2, 0);\r
+ }\r
+\r
+ const unsigned uTripleArrayBytes = TRIPLE_COUNT*sizeof(TripleCount);\r
+ TripleCounts = (TripleCount *) malloc(uTripleArrayBytes);\r
+ if (0 == TripleCounts)\r
+ Quit("Not enough memory (TripleCounts)");\r
+ memset(TripleCounts, 0, uTripleArrayBytes);\r
+\r
+ for (unsigned uWord = 0; uWord < TRIPLE_COUNT; ++uWord)\r
+ {\r
+ TripleCount &tc = *(TripleCounts + uWord);\r
+ const unsigned uBytes = uSeqCount*sizeof(short);\r
+ tc.m_Counts = (unsigned short *) malloc(uBytes);\r
+ memset(tc.m_Counts, 0, uBytes);\r
+ }\r
+\r
+ for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
+ {\r
+ Seq &s = *(v[uSeqIndex]);\r
+ const unsigned uSeqLength = s.Length();\r
+ for (unsigned uPos = 0; uPos < uSeqLength - 2; ++uPos)\r
+ {\r
+ const unsigned uLetter1 = CharToLetterEx(s[uPos]);\r
+ if (uLetter1 >= 20)\r
+ continue;\r
+ const unsigned uLetter2 = CharToLetterEx(s[uPos+1]);\r
+ if (uLetter2 >= 20)\r
+ continue;\r
+ const unsigned uLetter3 = CharToLetterEx(s[uPos+2]);\r
+ if (uLetter3 >= 20)\r
+ continue;\r
+\r
+ const unsigned uWord = uLetter1 + uLetter2*20 + uLetter3*20*20;\r
+ assert(uWord < TRIPLE_COUNT);\r
+\r
+ TripleCount &tc = *(TripleCounts + uWord);\r
+ const unsigned uOldCount = tc.m_Counts[uSeqIndex];\r
+ if (0 == uOldCount)\r
+ ++(tc.m_uSeqCount);\r
+\r
+ ++(tc.m_Counts[uSeqIndex]);\r
+ }\r
+ }\r
+\r
+#if TRACE\r
+ {\r
+ Log("TripleCounts\n");\r
+ unsigned uGrandTotal = 0;\r
+ for (unsigned uWord = 0; uWord < TRIPLE_COUNT; ++uWord)\r
+ {\r
+ const TripleCount &tc = *(TripleCounts + uWord);\r
+ if (0 == tc.m_uSeqCount)\r
+ continue;\r
+\r
+ const unsigned uLetter3 = uWord/(20*20);\r
+ const unsigned uLetter2 = (uWord - uLetter3*20*20)/20;\r
+ const unsigned uLetter1 = uWord%20;\r
+ Log("Word %6u %c%c%c %6u",\r
+ uWord,\r
+ LetterToCharAmino(uLetter1),\r
+ LetterToCharAmino(uLetter2),\r
+ LetterToCharAmino(uLetter3),\r
+ tc.m_uSeqCount);\r
+\r
+ unsigned uSeqCountWithThisWord = 0;\r
+ for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
+ {\r
+ const unsigned uCount = tc.m_Counts[uSeqIndex];\r
+ if (uCount > 0)\r
+ {\r
+ ++uSeqCountWithThisWord;\r
+ Log(" %u=%u", uSeqIndex, uCount);\r
+ uGrandTotal += uCount;\r
+ }\r
+ }\r
+ if (uSeqCountWithThisWord != tc.m_uSeqCount)\r
+ Log(" *** SQ ERROR *** %u %u", tc.m_uSeqCount, uSeqCountWithThisWord);\r
+ Log("\n");\r
+ }\r
+ \r
+ unsigned uTotalBySeqLength = 0;\r
+ for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
+ {\r
+ Seq &s = *(v[uSeqIndex]);\r
+ const unsigned uSeqLength = s.Length();\r
+ uTotalBySeqLength += uSeqLength - 2;\r
+ }\r
+ if (uGrandTotal != uTotalBySeqLength)\r
+ Log("*** TOTALS DISAGREE *** %u %u\n", uGrandTotal, uTotalBySeqLength);\r
+ }\r
+#endif\r
+\r
+ const unsigned uSeqListBytes = uSeqCount*sizeof(unsigned);\r
+ unsigned short *SeqList = (unsigned short *) malloc(uSeqListBytes);\r
+\r
+ for (unsigned uWord = 0; uWord < TRIPLE_COUNT; ++uWord)\r
+ {\r
+ const TripleCount &tc = *(TripleCounts + uWord);\r
+ if (0 == tc.m_uSeqCount)\r
+ continue;\r
+\r
+ unsigned uSeqCountFound = 0;\r
+ memset(SeqList, 0, uSeqListBytes);\r
+\r
+ for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
+ {\r
+ if (tc.m_Counts[uSeqIndex] > 0)\r
+ {\r
+ SeqList[uSeqCountFound] = uSeqIndex;\r
+ ++uSeqCountFound;\r
+ if (uSeqCountFound == tc.m_uSeqCount)\r
+ break;\r
+ }\r
+ }\r
+ assert(uSeqCountFound == tc.m_uSeqCount);\r
+\r
+ for (unsigned uSeq1 = 0; uSeq1 < uSeqCountFound; ++uSeq1)\r
+ {\r
+ const unsigned uSeqIndex1 = SeqList[uSeq1];\r
+ const unsigned uCount1 = tc.m_Counts[uSeqIndex1];\r
+ for (unsigned uSeq2 = 0; uSeq2 < uSeq1; ++uSeq2)\r
+ {\r
+ const unsigned uSeqIndex2 = SeqList[uSeq2];\r
+ const unsigned uCount2 = tc.m_Counts[uSeqIndex2];\r
+ const unsigned uMinCount = uCount1 < uCount2 ? uCount1 : uCount2;\r
+ const double d = DF.GetDist(uSeqIndex1, uSeqIndex2);\r
+ DF.SetDist(uSeqIndex1, uSeqIndex2, (float) (d + uMinCount));\r
+ }\r
+ }\r
+ }\r
+ delete[] SeqList;\r
+ free(TripleCounts);\r
+\r
+ unsigned uDone = 0;\r
+ const unsigned uTotal = (uSeqCount*(uSeqCount - 1))/2;\r
+ for (unsigned uSeq1 = 0; uSeq1 < uSeqCount; ++uSeq1)\r
+ {\r
+ DF.SetDist(uSeq1, uSeq1, 0.0);\r
+\r
+ const Seq &s1 = *(v[uSeq1]);\r
+ const unsigned uLength1 = s1.Length();\r
+\r
+ for (unsigned uSeq2 = 0; uSeq2 < uSeq1; ++uSeq2)\r
+ {\r
+ const Seq &s2 = *(v[uSeq2]);\r
+ const unsigned uLength2 = s2.Length();\r
+ unsigned uMinLength = uLength1 < uLength2 ? uLength1 : uLength2;\r
+ if (uMinLength < 3)\r
+ {\r
+ DF.SetDist(uSeq1, uSeq2, 1.0);\r
+ continue;\r
+ }\r
+\r
+ const double dTripleCount = DF.GetDist(uSeq1, uSeq2);\r
+ if (dTripleCount == 0)\r
+ {\r
+ DF.SetDist(uSeq1, uSeq2, 1.0);\r
+ continue;\r
+ }\r
+ double dNormalizedTripletScore = dTripleCount/(uMinLength - 2);\r
+ //double dEstimatedPairwiseIdentity = exp(0.3912*log(dNormalizedTripletScore));\r
+ //if (dEstimatedPairwiseIdentity > 1)\r
+ // dEstimatedPairwiseIdentity = 1;\r
+// DF.SetDist(uSeq1, uSeq2, (float) (1.0 - dEstimatedPairwiseIdentity));\r
+ DF.SetDist(uSeq1, uSeq2, (float) dNormalizedTripletScore);\r
+\r
+#if TRACE\r
+ {\r
+ Log("%s - %s Triplet count = %g Lengths %u, %u Estimated pwid = %g\n",\r
+ s1.GetName(), s2.GetName(), dTripleCount, uLength1, uLength2,\r
+ dEstimatedPairwiseIdentity);\r
+ }\r
+#endif\r
+ if (uDone%1000 == 0)\r
+ Progress(uDone, uTotal);\r
+ }\r
+ }\r
+ ProgressStepsDone();\r
+ }\r
git://source.jalview.org / jabaws.git / blobdiff