uchime_src/make3way.cpp

   1 #include "myutils.h"\r
   2 #include "sfasta.h"\r
   3 #include "path.h"\r
   4 #include "dp.h"\r
   5 \r
   6 void Make3Way(const SeqData &QSD, const SeqData &ASD, const SeqData &BSD,\r
   7   const string &PathQA, const string &PathQB,\r
   8   string &Q3, string &A3, string &B3)\r
   9         {\r
  10         Q3.clear();\r
  11         A3.clear();\r
  12         B3.clear();\r
  13 \r
  14 #if     DEBUG\r
  15         {\r
  16         unsigned QLen = 0;\r
  17         unsigned ALen = 0;\r
  18         for (unsigned i = 0; i < SIZE(PathQA); ++i)\r
  19                 {\r
  20                 char c = PathQA[i];\r
  21                 if (c == 'M' || c == 'D')\r
  22                         ++QLen;\r
  23                 if (c == 'M' || c == 'I')\r
  24                         ++ALen;\r
  25                 }\r
  26         asserta(QLen == QSD.L);\r
  27         asserta(ALen == ASD.L);\r
  28         }\r
  29         {\r
  30         unsigned QLen = 0;\r
  31         unsigned BLen = 0;\r
  32         for (unsigned i = 0; i < SIZE(PathQB); ++i)\r
  33                 {\r
  34                 char c = PathQB[i];\r
  35                 if (c == 'M' || c == 'D')\r
  36                         ++QLen;\r
  37                 if (c == 'M' || c == 'I')\r
  38                         ++BLen;\r
  39                 }\r
  40         asserta(QLen == QSD.L);\r
  41         asserta(BLen == BSD.L);\r
  42         }\r
  43 #endif\r
  44 \r
  45         const byte *Q = QSD.Seq;\r
  46         const byte *A = ASD.Seq;\r
  47         const byte *B = BSD.Seq;\r
  48 \r
  49         unsigned LQ = QSD.L;\r
  50         unsigned LA = ASD.L;\r
  51         unsigned LB = BSD.L;\r
  52 \r
  53         vector<unsigned> InsertCountsA(LQ+1, 0);\r
  54         unsigned QPos = 0;\r
  55         for (unsigned i = 0; i < SIZE(PathQA); ++i)\r
  56                 {\r
  57                 char c = PathQA[i];\r
  58                 if (c == 'M' || c == 'D')\r
  59                         ++QPos;\r
  60                 else\r
  61                         {\r
  62                         asserta(c == 'I');\r
  63                         asserta(QPos <= LQ);\r
  64                         ++(InsertCountsA[QPos]);\r
  65                         }\r
  66                 }\r
  67 \r
  68         vector<unsigned> InsertCountsB(LQ+1, 0);\r
  69         QPos = 0;\r
  70         for (unsigned i = 0; i < SIZE(PathQB); ++i)\r
  71                 {\r
  72                 char c = PathQB[i];\r
  73                 if (c == 'M' || c == 'D')\r
  74                         ++QPos;\r
  75                 else\r
  76                         {\r
  77                         asserta(c == 'I');\r
  78                         asserta(QPos <= LQ);\r
  79                         ++(InsertCountsB[QPos]);\r
  80                         }\r
  81                 }\r
  82 \r
  83         vector<unsigned> InsertCounts;\r
  84         for (unsigned i = 0; i <= LQ; ++i)\r
  85                 {\r
  86                 unsigned is = max(InsertCountsA[i], InsertCountsB[i]);\r
  87                 InsertCounts.push_back(is);\r
  88                 }\r
  89 \r
  90         for (unsigned i = 0; i < LQ; ++i)\r
  91                 {\r
  92                 for (unsigned k = 0; k < InsertCounts[i]; ++k)\r
  93                         Q3.push_back('-');\r
  94                 asserta(i < LQ);\r
  95                 Q3.push_back(toupper(Q[i]));\r
  96                 }\r
  97         for (unsigned k = 0; k < InsertCounts[LQ]; ++k)\r
  98                 Q3.push_back('-');\r
  99 \r
 100 // A\r
 101         QPos = 0;\r
 102         unsigned APos = 0;\r
 103         unsigned is = 0;\r
 104         for (unsigned i = 0; i < SIZE(PathQA); ++i)\r
 105                 {\r
 106                 char c = PathQA[i];\r
 107                 if (c == 'M' || c == 'D')\r
 108                         {\r
 109                         unsigned isq = InsertCounts[QPos];\r
 110                         asserta(is <= isq);\r
 111                         for (unsigned i = 0; i < InsertCounts[QPos]-is; ++i)\r
 112                                 A3.push_back('-');\r
 113                         is = 0;\r
 114                         ++QPos;\r
 115                         }\r
 116                 if (c == 'M')\r
 117                         {\r
 118                         asserta(APos < LA);\r
 119                         A3.push_back(toupper(A[APos++]));\r
 120                         }\r
 121                 else if (c == 'D')\r
 122                         A3.push_back('-');\r
 123                 else if (c == 'I')\r
 124                         {\r
 125                         ++is;\r
 126                         asserta(APos < LA);\r
 127                         A3.push_back(toupper(A[APos++]));\r
 128                         }\r
 129                 }\r
 130         asserta(is <= InsertCounts[LQ]);\r
 131         for (unsigned k = 0; k < InsertCounts[LQ]-is; ++k)\r
 132                 A3.push_back('-');\r
 133         asserta(QPos == LQ);\r
 134         asserta(APos == LA);\r
 135 \r
 136 // B\r
 137         QPos = 0;\r
 138         unsigned BPos = 0;\r
 139         is = 0;\r
 140         for (unsigned i = 0; i < SIZE(PathQB); ++i)\r
 141                 {\r
 142                 char c = PathQB[i];\r
 143                 if (c == 'M' || c == 'D')\r
 144                         {\r
 145                         asserta(is <= InsertCounts[QPos]);\r
 146                         for (unsigned i = 0; i < InsertCounts[QPos]-is; ++i)\r
 147                                 B3.push_back('-');\r
 148                         is = 0;\r
 149                         ++QPos;\r
 150                         }\r
 151                 if (c == 'M')\r
 152                         {\r
 153                         asserta(BPos < LB);\r
 154                         B3.push_back(toupper(B[BPos++]));\r
 155                         }\r
 156                 else if (c == 'D')\r
 157                         B3.push_back('-');\r
 158                 else if (c == 'I')\r
 159                         {\r
 160                         ++is;\r
 161                         asserta(BPos < LB);\r
 162                         B3.push_back(toupper(B[BPos++]));\r
 163                         }\r
 164                 }\r
 165         asserta(is <= InsertCounts[LQ]);\r
 166         for (unsigned k = 0; k < InsertCounts[LQ]-is; ++k)\r
 167                 B3.push_back('-');\r
 168         asserta(APos == LA);\r
 169         asserta(BPos == LB);\r
 170 \r
 171         asserta(SIZE(Q3) == SIZE(A3));\r
 172         asserta(SIZE(Q3) == SIZE(B3));\r
 173         }\r