uchime_src/alignchimel.cpp

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