alignchimel.cpp

   1 //uchime by Robert C. Edgar http://drive5.com/uchime This code is donated to the public domain.\r
   2 \r
   3 #include "myutils.h"\r
   4 #include "seq.h"\r
   5 #include "chime.h"\r
   6 \r
   7 #define TRACE   0\r
   8 \r
   9 /***\r
  10 Let:\r
  11         S[i] =  Score of col i: 0=no SNP, +1 = Y, -3 = N or A.\r
  12 \r
  13         V[k] =  Best segment score from j, j+1 .. k for all possible j\r
  14                         max(j) Sum i=j..k S[i]\r
  15 \r
  16 Recursion relation:\r
  17         V[k] =  S[k] + max (V[k-1], 0)\r
  18 ***/\r
  19 \r
  20 void AlignChimeGlobal3(const string &Q3, const string &A3, const string &B3,\r
  21   const string &QLabel, const string &ALabel, const string &BLabel,\r
  22   ChimeHit2 &Hit);\r
  23 \r
  24 void Make3Way(const SeqData &SDQ, const SeqData &SDA, const SeqData &SDB,\r
  25   const string &PathQA, const string &PathQB,\r
  26   string &Q3, string &A3, string &B3);\r
  27 \r
  28 double GetScore2(double Y, double N, double A);\r
  29 \r
  30 void AlignChimeLocal3(const string &Q3, const string &A3, const string &B3,\r
  31   const string &QLabel, const string &ALabel, const string &BLabel,\r
  32   ChimeHit2 &Hit)\r
  33         {\r
  34         Hit.Clear();\r
  35 \r
  36         const byte *Q3Seq = (const byte *) Q3.c_str();\r
  37         const byte *A3Seq = (const byte *) A3.c_str();\r
  38         const byte *B3Seq = (const byte *) B3.c_str();\r
  39 \r
  40         const unsigned ColCount = SIZE(Q3);\r
  41         asserta(SIZE(A3) == ColCount && SIZE(B3) == ColCount);\r
  42 \r
  43         vector<float> ColScoresA(ColCount, 0.0f);\r
  44         vector<float> ColScoresB(ColCount, 0.0f);\r
  45 \r
  46         float ScoreN = -(float) opt_xn;\r
  47         unsigned QL = 0;\r
  48         for (unsigned Col = 0; Col < ColCount; ++Col)\r
  49                 {\r
  50                 char q = Q3Seq[Col];\r
  51                 char a = A3Seq[Col];\r
  52                 char b = B3Seq[Col];\r
  53 \r
  54                 if (!isgap(q))\r
  55                         ++QL;\r
  56 \r
  57                 if (q == a && q == b && a == b)\r
  58                         continue;\r
  59 \r
  60                 if (isgap(q) || isgap(a) || isgap(b))\r
  61                         continue;\r
  62 \r
  63                 if (Col > 0 && (isgap(Q3Seq[Col-1]) || isgap(A3Seq[Col-1]) || isgap(B3Seq[Col-1])))\r
  64                         continue;\r
  65 \r
  66                 if (Col + 1 < ColCount && (isgap(Q3Seq[Col+1]) || isgap(A3Seq[Col+1]) || isgap(B3Seq[Col+1])))\r
  67                         continue;\r
  68 \r
  69                 if (q == a && q != b)\r
  70                         ColScoresA[Col] = 1;\r
  71                 else\r
  72                         ColScoresA[Col] = ScoreN;\r
  73 \r
  74                 if (q == b && q != a)\r
  75                         ColScoresB[Col] = 1;\r
  76                 else\r
  77                         ColScoresB[Col] = ScoreN;\r
  78                 }\r
  79 \r
  80         vector<float> LVA(ColCount, 0.0f);\r
  81         vector<float> LVB(ColCount, 0.0f);\r
  82 \r
  83         LVA[0] = ColScoresA[0];\r
  84         LVB[0] = ColScoresB[0];\r
  85         for (unsigned Col = 1; Col < ColCount; ++Col)\r
  86                 {\r
  87                 LVA[Col] = max(LVA[Col-1], 0.0f) + ColScoresA[Col];\r
  88                 LVB[Col] = max(LVB[Col-1], 0.0f) + ColScoresB[Col];\r
  89                 }\r
  90 \r
  91         vector<float> RVA(ColCount, 0.0f);\r
  92         vector<float> RVB(ColCount, 0.0f);\r
  93 \r
  94         RVA[ColCount-1] = ColScoresA[ColCount-1];\r
  95         RVB[ColCount-1] = ColScoresB[ColCount-1];\r
  96         for (int Col = ColCount-2; Col >= 0; --Col)\r
  97                 {\r
  98                 RVA[Col] = max(RVA[Col+1], 0.0f) + ColScoresA[Col];\r
  99                 RVB[Col] = max(RVB[Col+1], 0.0f) + ColScoresB[Col];\r
 100                 }\r
 101 \r
 102         bool FirstA = true;\r
 103         float MaxSum = 0.0;\r
 104         unsigned ColX = UINT_MAX;\r
 105         for (unsigned Col = 1; Col < ColCount-1; ++Col)\r
 106                 {\r
 107                 float Sum = LVA[Col] + RVB[Col+1];\r
 108                 if (Sum > MaxSum)\r
 109                         {\r
 110                         FirstA = true;\r
 111                         MaxSum = Sum;\r
 112                         ColX = Col;\r
 113                         }\r
 114                 }\r
 115 \r
 116         for (unsigned Col = 1; Col < ColCount-1; ++Col)\r
 117                 {\r
 118                 float Sum = LVB[Col] + RVA[Col+1];\r
 119                 if (Sum > MaxSum)\r
 120                         {\r
 121                         FirstA = false;\r
 122                         MaxSum = Sum;\r
 123                         ColX = Col;\r
 124                         }\r
 125                 }\r
 126         if (ColX == UINT_MAX)\r
 127                 return;\r
 128 \r
 129         unsigned ColLo = UINT_MAX;\r
 130         unsigned ColHi = UINT_MAX;\r
 131         if (FirstA)\r
 132                 {\r
 133                 float Sum = 0.0f;\r
 134                 for (int Col = ColX; Col >= 0; --Col)\r
 135                         {\r
 136                         Sum += ColScoresA[Col];\r
 137                         if (Sum >= LVA[ColX])\r
 138                                 {\r
 139                                 ColLo = Col;\r
 140                                 break;\r
 141                                 }\r
 142                         }\r
 143                 asserta(Sum >= LVA[ColX]);\r
 144                 Sum = 0.0f;\r
 145                 for (unsigned Col = ColX+1; Col < ColCount; ++Col)\r
 146                         {\r
 147                         Sum += ColScoresB[Col];\r
 148                         if (Sum >= RVB[ColX])\r
 149                                 {\r
 150                                 ColHi = Col;\r
 151                                 break;\r
 152                                 }\r
 153                         }\r
 154                 asserta(Sum >= RVB[ColX]);\r
 155                 }\r
 156         else\r
 157                 {\r
 158                 float Sum = 0.0f;\r
 159                 for (int Col = ColX; Col >= 0; --Col)\r
 160                         {\r
 161                         Sum += ColScoresB[Col];\r
 162                         if (Sum >= LVB[ColX])\r
 163                                 {\r
 164                                 ColLo = Col;\r
 165                                 break;\r
 166                                 }\r
 167                         }\r
 168                 asserta(Sum >= LVB[ColX]);\r
 169                 Sum = 0.0f;\r
 170                 for (unsigned Col = ColX+1; Col < ColCount; ++Col)\r
 171                         {\r
 172                         Sum += ColScoresA[Col];\r
 173                         if (Sum >= RVA[ColX])\r
 174                                 {\r
 175                                 ColHi = Col;\r
 176                                 break;\r
 177                                 }\r
 178                         }\r
 179                 asserta(Sum >= RVA[ColX]);\r
 180                 }\r
 181 \r
 182         unsigned ColXHi = ColX;\r
 183         for (unsigned Col = ColX + 1; Col < ColCount; ++Col)\r
 184                 {\r
 185                 char q = Q3Seq[Col];\r
 186                 char a = A3Seq[Col];\r
 187                 char b = B3Seq[Col];\r
 188                 \r
 189                 if (q == a && q == b && !isgap(q))\r
 190                         ColXHi = Col;\r
 191                 else\r
 192                         break;\r
 193                 }\r
 194 \r
 195         unsigned ColXLo = ColX;\r
 196         for (int Col = (int) ColX - 1; Col >= 0; --Col)\r
 197                 {\r
 198                 char q = Q3Seq[Col];\r
 199                 char a = A3Seq[Col];\r
 200                 char b = B3Seq[Col];\r
 201                 \r
 202                 if (q == a && q == b && !isgap(q))\r
 203                         ColXLo = Col;\r
 204                 else\r
 205                         break;\r
 206                 }\r
 207 \r
 208         unsigned IdQA = 0;\r
 209         unsigned IdQB = 0;\r
 210         unsigned IdAB = 0;\r
 211         unsigned NQA = 0;\r
 212         unsigned NQB = 0;\r
 213         unsigned NAB = 0;\r
 214         for (unsigned Col = 0; Col < ColCount; ++Col)\r
 215                 {\r
 216                 char q = Q3Seq[Col];\r
 217                 char a = A3Seq[Col];\r
 218                 char b = B3Seq[Col];\r
 219 \r
 220                 if (!isgap(q) && !isgap(a))\r
 221                         {\r
 222                         ++NQA;\r
 223                         if (q == a)\r
 224                                 ++IdQA;\r
 225                         }\r
 226 \r
 227                 if (!isgap(q) && !isgap(b))\r
 228                         {\r
 229                         ++NQB;\r
 230                         if (q == b)\r
 231                                 ++IdQB;\r
 232                         }\r
 233 \r
 234                 if (!isgap(a) && !isgap(b))\r
 235                         {\r
 236                         ++NAB;\r
 237                         if (a == b)\r
 238                                 ++IdAB;\r
 239                         }\r
 240                 }\r
 241 \r
 242         Hit.PctIdQA = Pct(IdQA, NQA);\r
 243         Hit.PctIdQB = Pct(IdQB, NQB);\r
 244         Hit.PctIdAB = Pct(IdAB, NAB);\r
 245 \r
 246         unsigned LIdQA = 0;\r
 247         unsigned LIdQB = 0;\r
 248         for (unsigned Col = ColLo; Col < ColXLo; ++Col)\r
 249                 {\r
 250                 char q = Q3Seq[Col];\r
 251                 char a = A3Seq[Col];\r
 252                 char b = B3Seq[Col];\r
 253 \r
 254                 if (!isgap(q) && !isgap(a))\r
 255                         {\r
 256                         if (q == a)\r
 257                                 ++LIdQA;\r
 258                         }\r
 259 \r
 260                 if (!isgap(q) && !isgap(b))\r
 261                         {\r
 262                         if (q == b)\r
 263                                 ++LIdQB;\r
 264                         }\r
 265                 }\r
 266 \r
 267         unsigned RIdQA = 0;\r
 268         unsigned RIdQB = 0;\r
 269         for (unsigned Col = ColXHi+1; Col <= ColHi; ++Col)\r
 270                 {\r
 271                 char q = Q3Seq[Col];\r
 272                 char a = A3Seq[Col];\r
 273                 char b = B3Seq[Col];\r
 274 \r
 275                 if (!isgap(q) && !isgap(a))\r
 276                         {\r
 277                         if (q == a)\r
 278                                 ++RIdQA;\r
 279                         }\r
 280 \r
 281                 if (!isgap(q) && !isgap(b))\r
 282                         {\r
 283                         if (q == b)\r
 284                                 ++RIdQB;\r
 285                         }\r
 286                 }\r
 287 \r
 288         unsigned IdDiffL = max(LIdQA, LIdQB) - min(LIdQA, LIdQB);\r
 289         unsigned IdDiffR = max(RIdQA, RIdQB) - min(RIdQA, RIdQB);\r
 290         unsigned MinIdDiff = min(IdDiffL, IdDiffR);\r
 291         unsigned ColRange = ColHi - ColLo + 1;\r
 292         if (opt_queryfract > 0.0f && float(ColRange)/float(QL) < opt_queryfract)\r
 293                 return;\r
 294 \r
 295 //      double Div = Pct(MinIdDiff, QSD.L);\r
 296 \r
 297 #if     TRACE\r
 298         {\r
 299         Log("  Col  A Q B   ScoreA   ScoreB      LVA      LVB      RVA      RVB\n");\r
 300         Log("-----  - - -  -------  -------  -------  -------  -------  -------\n");\r
 301         for (unsigned Col = 0; Col < ColCount; ++Col)\r
 302                 {\r
 303                 if (ColScoresA[Col] == 0.0 && ColScoresB[Col] == 0.0)\r
 304                         continue;\r
 305 \r
 306                 char q = Q3Seq[Col];\r
 307                 char a = A3Seq[Col];\r
 308                 char b = B3Seq[Col];\r
 309                 Log("%5u  %c %c %c", Col, a, q, b);\r
 310 \r
 311                 if (ColScoresA[Col] == 0.0)\r
 312                         Log("  %7.7s", "");\r
 313                 else\r
 314                         Log("  %7.1f", ColScoresA[Col]);\r
 315 \r
 316                 if (ColScoresB[Col] == 0.0)\r
 317                         Log("  %7.7s", "");\r
 318                 else\r
 319                         Log("  %7.1f", ColScoresB[Col]);\r
 320 \r
 321                 Log("  %7.1f  %7.1f  %7.1f  %7.1f", LVA[Col], LVB[Col], RVA[Col], RVB[Col]);\r
 322 \r
 323                 Log("\n");\r
 324                 }\r
 325         Log("\n");\r
 326         Log("MaxSum %.1f, ColLo %u, ColXLo %u, ColX %u, ColXHi %u, ColHi %u, AF %c\n",\r
 327           MaxSum, ColLo, ColXLo, ColX, ColXHi, ColHi, tof(FirstA));\r
 328         Log("  LIdQA %u, LIdQB %u, RIdQA %u, RIdQB %u\n", LIdQA, LIdQB, RIdQA, RIdQB);\r
 329         }\r
 330 #endif\r
 331 \r
 332         string Q3L;\r
 333         string A3L;\r
 334         string B3L;\r
 335         for (unsigned Col = ColLo; Col <= ColHi; ++Col)\r
 336                 {\r
 337                 char q = Q3[Col];\r
 338                 char a = A3[Col];\r
 339                 char b = B3[Col];\r
 340 \r
 341                 Q3L += q;\r
 342                 A3L += a;\r
 343                 B3L += b;\r
 344                 }\r
 345 \r
 346         AlignChimeGlobal3(Q3L, A3L, B3L, QLabel, ALabel, BLabel, Hit);\r
 347 \r
 348 #if     0\r
 349 // CS SNPs\r
 350         Hit.CS_LY = 0;\r
 351         Hit.CS_LN = 0;\r
 352         Hit.CS_RY = 0;\r
 353         Hit.CS_RN = 0;\r
 354         Hit.CS_LA = 0;\r
 355         Hit.CS_RA = 0;\r
 356         for (unsigned Col = ColLo; Col <= ColHi; ++Col)\r
 357                 {\r
 358                 char q = Q3Seq[Col];\r
 359                 char a = A3Seq[Col];\r
 360                 char b = B3Seq[Col];\r
 361                 if (q == a && q == b && a == b)\r
 362                         continue;\r
 363                 if (isgap(q) || isgap(a) || isgap(b))\r
 364                         continue;\r
 365                 if (Col > 0 && (isgap(Q3Seq[Col-1]) || isgap(A3Seq[Col-1]) || isgap(B3Seq[Col-1])))\r
 366                         continue;\r
 367                 if (Col + 1 < ColCount && (isgap(Q3Seq[Col+1]) || isgap(A3Seq[Col+1]) || isgap(B3Seq[Col+1])))\r
 368                         continue;\r
 369 \r
 370                 if (!FirstA)\r
 371                         swap(a, b);\r
 372 \r
 373                 if (Col < ColXLo)\r
 374                         {\r
 375                         if (q == a && q != b)\r
 376                                 ++Hit.CS_LY;\r
 377                         else if (q == b && q != a)\r
 378                                 ++Hit.CS_LN;\r
 379                         else\r
 380                                 ++Hit.CS_LA;\r
 381                         }\r
 382                 else if (Col > ColXHi)\r
 383                         {\r
 384                         if (q == b && q != a)\r
 385                                 ++Hit.CS_RY;\r
 386                         else if (q == a && q != b)\r
 387                                 ++Hit.CS_RN;\r
 388                         else\r
 389                                 ++Hit.CS_RA;\r
 390                         }\r
 391                 }\r
 392 \r
 393         double ScoreL = GetScore2(Hit.CS_LY, Hit.CS_LN, Hit.CS_LA);\r
 394         double ScoreR = GetScore2(Hit.CS_RY, Hit.CS_RN, Hit.CS_RA);\r
 395         Hit.Score = ScoreL*ScoreR;\r
 396 \r
 397         //Hit.QSD = QSD;\r
 398         //if (FirstA)\r
 399         //      {\r
 400         //      Hit.ASD = ASD;\r
 401         //      Hit.BSD = BSD;\r
 402         //      Hit.PathQA = PathQA;\r
 403         //      Hit.PathQB = PathQB;\r
 404         //      }\r
 405         //else\r
 406         //      {\r
 407         //      Hit.ASD = BSD;\r
 408         //      Hit.BSD = ASD;\r
 409         //      }\r
 410 \r
 411         //Hit.ColLo = ColLo;\r
 412         //Hit.ColXLo = ColXLo;\r
 413         //Hit.ColXHi = ColXHi;\r
 414         //Hit.ColHi = ColHi;\r
 415         //Hit.Div = Div;\r
 416 \r
 417 //      Hit.LogMe();\r
 418 #endif\r
 419         }\r