--- /dev/null
+#include "myutils.h"\r
+#include "ultra.h"\r
+#include "chime.h"\r
+#include "uc.h"\r
+#include "dp.h"\r
+#include <set>\r
+#include <algorithm>\r
+\r
+#define TRACE 0\r
+\r
+extern FILE *g_fUChime;\r
+\r
+void GetCandidateParents(Ultra &U, const SeqData &QSD, float AbQ,\r
+ vector<unsigned> &Parents);\r
+\r
+void AlignChime(const SeqData &QSD, const SeqData &ASD, const SeqData &BSD,\r
+ const string &PathQA, const string &PathQB, ChimeHit2 &Hit);\r
+\r
+double GetFractIdGivenPath(const byte *A, const byte *B, const char *Path, bool Nucleo);\r
+\r
+static void GetSmoothedIdVec(const SeqData &QSD, const SeqData &PSD, const string &Path,\r
+ vector<unsigned> &IdVec, unsigned d)\r
+ {\r
+ IdVec.clear();\r
+ const unsigned ColCount = SIZE(Path);\r
+\r
+ const byte *Q = QSD.Seq;\r
+ const byte *P = PSD.Seq;\r
+\r
+ const unsigned QL = QSD.L;\r
+ const unsigned PL = PSD.L;\r
+\r
+ if (QL <= d)\r
+ {\r
+ IdVec.resize(QSD.L, 0);\r
+ return;\r
+ }\r
+\r
+ unsigned QPos = 0;\r
+ unsigned PPos = 0;\r
+\r
+ vector<bool> SameVec;\r
+ SameVec.reserve(QL);\r
+ for (unsigned Col = 0; Col < ColCount; ++Col)\r
+ {\r
+ char c = Path[Col];\r
+\r
+ bool Same = false;\r
+ if (c == 'M')\r
+ {\r
+ byte q = Q[QPos];\r
+ byte p = P[PPos];\r
+ Same = (toupper(q) == toupper(p));\r
+ }\r
+\r
+ if (c == 'M' || c == 'D')\r
+ {\r
+ ++QPos;\r
+ SameVec.push_back(Same);\r
+ }\r
+\r
+ if (c == 'M' || c == 'I')\r
+ ++PPos;\r
+ }\r
+\r
+ asserta(SIZE(SameVec) == QL);\r
+\r
+ unsigned n = 0;\r
+ for (unsigned QPos = 0; QPos < d; ++QPos)\r
+ {\r
+ if (SameVec[QPos])\r
+ ++n;\r
+ IdVec.push_back(n);\r
+ }\r
+\r
+ for (unsigned QPos = d; QPos < QL; ++QPos)\r
+ {\r
+ if (SameVec[QPos])\r
+ ++n;\r
+ IdVec.push_back(n);\r
+ if (SameVec[QPos-d])\r
+ --n;\r
+ }\r
+ asserta(SIZE(IdVec) == QL);\r
+\r
+#if TRACE\r
+ {\r
+ Log("\n");\r
+ Log("GetSmoothedIdVec\n");\r
+ unsigned QPos = 0;\r
+ unsigned PPos = 0;\r
+ Log("Q P Same Id\n");\r
+ Log("- - ---- -------\n");\r
+ for (unsigned Col = 0; Col < ColCount; ++Col)\r
+ {\r
+ char c = Path[Col];\r
+\r
+ bool Same = false;\r
+ if (c == 'M')\r
+ {\r
+ byte q = Q[QPos];\r
+ byte p = P[PPos];\r
+ Same = (toupper(q) == toupper(p));\r
+ Log("%c %c %4c %7d\n", q, p, tof(Same), IdVec[QPos]);\r
+ }\r
+\r
+ if (c == 'M' || c == 'D')\r
+ ++QPos;\r
+ if (c == 'M' || c == 'I')\r
+ ++PPos;\r
+ }\r
+ }\r
+#endif\r
+ }\r
+\r
+bool SearchChime(Ultra &U, const SeqData &QSD, float QAb, \r
+ const AlnParams &AP, const AlnHeuristics &AH, HSPFinder &HF,\r
+ float MinFractId, ChimeHit2 &Hit)\r
+ {\r
+ Hit.Clear();\r
+ Hit.QLabel = QSD.Label;\r
+\r
+ if (opt_verbose)\r
+ {\r
+ Log("\n");\r
+ Log("SearchChime()\n");\r
+ Log("Query>%s\n", QSD.Label);\r
+ }\r
+\r
+ vector<unsigned> Parents;\r
+ GetCandidateParents(U, QSD, QAb, Parents);\r
+\r
+ unsigned ParentCount = SIZE(Parents);\r
+ if (ParentCount <= 1)\r
+ {\r
+ if (opt_verbose)\r
+ Log("%u candidate parents, done.\n", ParentCount);\r
+ return false;\r
+ }\r
+\r
+ if (opt_fastalign)\r
+ HF.SetA(QSD);\r
+ HSPFinder *ptrHF = (opt_fastalign ? &HF : 0);\r
+\r
+ unsigned ChunkLength;\r
+ vector<unsigned> ChunkLos;\r
+ GetChunkInfo(QSD.L, ChunkLength, ChunkLos);\r
+ const unsigned ChunkCount = SIZE(ChunkLos);\r
+\r
+ vector<unsigned> ChunkIndexToBestId(ChunkCount, 0);\r
+ vector<unsigned> ChunkIndexToBestParentIndex(ChunkCount, UINT_MAX);\r
+\r
+ vector<SeqData> PSDs;\r
+ vector<string> Paths;\r
+ double TopPctId = 0.0;\r
+ unsigned TopParentIndex = UINT_MAX;\r
+ unsigned QL = QSD.L;\r
+ vector<unsigned> MaxIdVec(QL, 0);\r
+ for (unsigned ParentIndex = 0; ParentIndex < ParentCount; ++ParentIndex)\r
+ {\r
+ unsigned ParentSeqIndex = Parents[ParentIndex];\r
+\r
+ SeqData PSD;\r
+ //PSD.Label = U.GetSeedLabel(ParentSeqIndex);\r
+ //PSD.Seq = U.GetSeedSeq(ParentSeqIndex);\r
+ //PSD.L = U.GetSeedLength(ParentSeqIndex);\r
+ //PSD.Index = ParentSeqIndex;\r
+ U.GetSeqData(ParentSeqIndex, PSD);\r
+ PSDs.push_back(PSD);\r
+\r
+ if (opt_fastalign)\r
+ HF.SetB(PSD);\r
+\r
+ PathData PD;\r
+\r
+ float HSPId;\r
+ bool Found = GlobalAlign(QSD, PSD, AP, AH, *ptrHF, MinFractId, HSPId, PD);\r
+ if (!Found)\r
+ {\r
+ Paths.push_back(""); \r
+ continue;\r
+ }\r
+\r
+ double PctId = 100.0*GetFractIdGivenPath(QSD.Seq, PSD.Seq, PD.Start, true);\r
+ if (opt_selfid && PctId == 100.0)\r
+ {\r
+ Paths.push_back(""); \r
+ continue;\r
+ }\r
+\r
+ if (PctId > TopPctId)\r
+ {\r
+ TopParentIndex = ParentIndex;\r
+ TopPctId = PctId;\r
+ if (TopPctId >= 100.0 - opt_mindiv)\r
+ {\r
+ if (opt_verbose)\r
+ {\r
+ Log(" %.1f%% >%s\n", TopPctId, PSD.Label);\r
+ Log(" Top hit exceeds ctl threshold, done.\n");\r
+ return false;\r
+ }\r
+ }\r
+ }\r
+\r
+ string Path = PD.Start;\r
+ Paths.push_back(Path);\r
+\r
+ vector<unsigned> IdVec;\r
+ GetSmoothedIdVec(QSD, PSD, Path, IdVec, opt_idsmoothwindow);\r
+\r
+ for (unsigned QPos = 0; QPos < QL; ++QPos)\r
+ if (IdVec[QPos] > MaxIdVec[QPos])\r
+ MaxIdVec[QPos] = IdVec[QPos];\r
+ }\r
+\r
+ vector<unsigned> BestParents;\r
+ for (unsigned k = 0; k < opt_maxp; ++k)\r
+ {\r
+ unsigned BestParent = UINT_MAX;\r
+ unsigned BestCov = 0;\r
+ for (unsigned ParentIndex = 0; ParentIndex < ParentCount; ++ParentIndex)\r
+ {\r
+ const SeqData &PSD = PSDs[ParentIndex];\r
+ const string &Path = Paths[ParentIndex];\r
+ if (Path == "")\r
+ continue;\r
+\r
+ vector<unsigned> IdVec;\r
+ GetSmoothedIdVec(QSD, PSD, Path, IdVec, opt_idsmoothwindow);\r
+\r
+ unsigned Cov = 0;\r
+ for (unsigned QPos = 0; QPos < QL; ++QPos)\r
+ if (IdVec[QPos] == MaxIdVec[QPos])\r
+ ++Cov;\r
+\r
+ if (Cov > BestCov)\r
+ {\r
+ BestParent = ParentIndex;\r
+ BestCov = Cov;\r
+ }\r
+ }\r
+\r
+ if (BestParent == UINT_MAX)\r
+ break;\r
+\r
+ BestParents.push_back(BestParent);\r
+ vector<unsigned> IdVec;\r
+\r
+ const SeqData &PSD = PSDs[BestParent];\r
+ const string &Path = Paths[BestParent];\r
+ GetSmoothedIdVec(QSD, PSD, Path, IdVec, opt_idsmoothwindow);\r
+ for (unsigned QPos = 0; QPos < QL; ++QPos)\r
+ if (IdVec[QPos] == MaxIdVec[QPos])\r
+ MaxIdVec[QPos] = UINT_MAX;\r
+ }\r
+\r
+ unsigned BestParentCount = SIZE(BestParents);\r
+\r
+ if (opt_verbose)\r
+ {\r
+ Log("%u/%u best parents\n", BestParentCount, ParentCount);\r
+ for (unsigned k = 0; k < BestParentCount; ++k)\r
+ {\r
+ unsigned i = BestParents[k];\r
+ Log(" %s\n", PSDs[i].Label);\r
+ }\r
+ }\r
+\r
+ bool Found = false;\r
+ for (unsigned k1 = 0; k1 < BestParentCount; ++k1)\r
+ {\r
+ unsigned i1 = BestParents[k1];\r
+ asserta(i1 < ParentCount);\r
+\r
+ const SeqData &PSD1 = PSDs[i1];\r
+ const string &Path1 = Paths[i1];\r
+\r
+ for (unsigned k2 = k1 + 1; k2 < BestParentCount; ++k2)\r
+ {\r
+ unsigned i2 = BestParents[k2];\r
+ asserta(i2 < ParentCount);\r
+ asserta(i2 != i1);\r
+\r
+ const SeqData &PSD2 = PSDs[i2];\r
+ const string &Path2 = Paths[i2];\r
+\r
+ ChimeHit2 Hit2;\r
+ AlignChime(QSD, PSD1, PSD2, Path1, Path2, Hit2);\r
+ Hit2.PctIdQT = TopPctId;\r
+\r
+ if (Hit2.Accept())\r
+ Found = true;\r
+\r
+ if (Hit2.Score > Hit.Score)\r
+ Hit = Hit2;\r
+\r
+ if (opt_verbose)\r
+ Hit2.LogMe();\r
+ }\r
+ }\r
+\r
+ return Found;\r
+ }\r
projects / mothur.git / blobdiff