added uchime_src folder. added biom parameter to make.shared. added biom as a current...

[mothur.git] / uchime_src / make3way.cpp

diff --git a/uchime_src/make3way.cpp b/uchime_src/make3way.cpp
new file mode 100644 (file)
index 0000000..ce88f86
--- /dev/null
+++ b/uchime_src/make3way.cpp
@@ -0,0 +1,173 @@
+#include "myutils.h"\r
+#include "sfasta.h"\r
+#include "path.h"\r
+#include "dp.h"\r
+\r
+void Make3Way(const SeqData &QSD, const SeqData &ASD, const SeqData &BSD,\r
+  const string &PathQA, const string &PathQB,\r
+  string &Q3, string &A3, string &B3)\r
+       {\r
+       Q3.clear();\r
+       A3.clear();\r
+       B3.clear();\r
+\r
+#if    DEBUG\r
+       {\r
+       unsigned QLen = 0;\r
+       unsigned ALen = 0;\r
+       for (unsigned i = 0; i < SIZE(PathQA); ++i)\r
+               {\r
+               char c = PathQA[i];\r
+               if (c == 'M' || c == 'D')\r
+                       ++QLen;\r
+               if (c == 'M' || c == 'I')\r
+                       ++ALen;\r
+               }\r
+       asserta(QLen == QSD.L);\r
+       asserta(ALen == ASD.L);\r
+       }\r
+       {\r
+       unsigned QLen = 0;\r
+       unsigned BLen = 0;\r
+       for (unsigned i = 0; i < SIZE(PathQB); ++i)\r
+               {\r
+               char c = PathQB[i];\r
+               if (c == 'M' || c == 'D')\r
+                       ++QLen;\r
+               if (c == 'M' || c == 'I')\r
+                       ++BLen;\r
+               }\r
+       asserta(QLen == QSD.L);\r
+       asserta(BLen == BSD.L);\r
+       }\r
+#endif\r
+\r
+       const byte *Q = QSD.Seq;\r
+       const byte *A = ASD.Seq;\r
+       const byte *B = BSD.Seq;\r
+\r
+       unsigned LQ = QSD.L;\r
+       unsigned LA = ASD.L;\r
+       unsigned LB = BSD.L;\r
+\r
+       vector<unsigned> InsertCountsA(LQ+1, 0);\r
+       unsigned QPos = 0;\r
+       for (unsigned i = 0; i < SIZE(PathQA); ++i)\r
+               {\r
+               char c = PathQA[i];\r
+               if (c == 'M' || c == 'D')\r
+                       ++QPos;\r
+               else\r
+                       {\r
+                       asserta(c == 'I');\r
+                       asserta(QPos <= LQ);\r
+                       ++(InsertCountsA[QPos]);\r
+                       }\r
+               }\r
+\r
+       vector<unsigned> InsertCountsB(LQ+1, 0);\r
+       QPos = 0;\r
+       for (unsigned i = 0; i < SIZE(PathQB); ++i)\r
+               {\r
+               char c = PathQB[i];\r
+               if (c == 'M' || c == 'D')\r
+                       ++QPos;\r
+               else\r
+                       {\r
+                       asserta(c == 'I');\r
+                       asserta(QPos <= LQ);\r
+                       ++(InsertCountsB[QPos]);\r
+                       }\r
+               }\r
+\r
+       vector<unsigned> InsertCounts;\r
+       for (unsigned i = 0; i <= LQ; ++i)\r
+               {\r
+               unsigned is = max(InsertCountsA[i], InsertCountsB[i]);\r
+               InsertCounts.push_back(is);\r
+               }\r
+\r
+       for (unsigned i = 0; i < LQ; ++i)\r
+               {\r
+               for (unsigned k = 0; k < InsertCounts[i]; ++k)\r
+                       Q3.push_back('-');\r
+               asserta(i < LQ);\r
+               Q3.push_back(toupper(Q[i]));\r
+               }\r
+       for (unsigned k = 0; k < InsertCounts[LQ]; ++k)\r
+               Q3.push_back('-');\r
+\r
+// A\r
+       QPos = 0;\r
+       unsigned APos = 0;\r
+       unsigned is = 0;\r
+       for (unsigned i = 0; i < SIZE(PathQA); ++i)\r
+               {\r
+               char c = PathQA[i];\r
+               if (c == 'M' || c == 'D')\r
+                       {\r
+                       unsigned isq = InsertCounts[QPos];\r
+                       asserta(is <= isq);\r
+                       for (unsigned i = 0; i < InsertCounts[QPos]-is; ++i)\r
+                               A3.push_back('-');\r
+                       is = 0;\r
+                       ++QPos;\r
+                       }\r
+               if (c == 'M')\r
+                       {\r
+                       asserta(APos < LA);\r
+                       A3.push_back(toupper(A[APos++]));\r
+                       }\r
+               else if (c == 'D')\r
+                       A3.push_back('-');\r
+               else if (c == 'I')\r
+                       {\r
+                       ++is;\r
+                       asserta(APos < LA);\r
+                       A3.push_back(toupper(A[APos++]));\r
+                       }\r
+               }\r
+       asserta(is <= InsertCounts[LQ]);\r
+       for (unsigned k = 0; k < InsertCounts[LQ]-is; ++k)\r
+               A3.push_back('-');\r
+       asserta(QPos == LQ);\r
+       asserta(APos == LA);\r
+\r
+// B\r
+       QPos = 0;\r
+       unsigned BPos = 0;\r
+       is = 0;\r
+       for (unsigned i = 0; i < SIZE(PathQB); ++i)\r
+               {\r
+               char c = PathQB[i];\r
+               if (c == 'M' || c == 'D')\r
+                       {\r
+                       asserta(is <= InsertCounts[QPos]);\r
+                       for (unsigned i = 0; i < InsertCounts[QPos]-is; ++i)\r
+                               B3.push_back('-');\r
+                       is = 0;\r
+                       ++QPos;\r
+                       }\r
+               if (c == 'M')\r
+                       {\r
+                       asserta(BPos < LB);\r
+                       B3.push_back(toupper(B[BPos++]));\r
+                       }\r
+               else if (c == 'D')\r
+                       B3.push_back('-');\r
+               else if (c == 'I')\r
+                       {\r
+                       ++is;\r
+                       asserta(BPos < LB);\r
+                       B3.push_back(toupper(B[BPos++]));\r
+                       }\r
+               }\r
+       asserta(is <= InsertCounts[LQ]);\r
+       for (unsigned k = 0; k < InsertCounts[LQ]-is; ++k)\r
+               B3.push_back('-');\r
+       asserta(APos == LA);\r
+       asserta(BPos == LB);\r
+\r
+       asserta(SIZE(Q3) == SIZE(A3));\r
+       asserta(SIZE(Q3) == SIZE(B3));\r
+       }\r