uchime_src/alignchime.cpp

   1 #include "myutils.h"\r
   2 #include "seq.h"\r
   3 #include "chime.h"\r
   4 #include "dp.h"\r
   5 \r
   6 #define TRACE           0\r
   7 #define TRACE_BS        0\r
   8 \r
   9 void Make3Way(const SeqData &SDQ, const SeqData &SDA, const SeqData &SDB,\r
  10   const string &PathQA, const string &PathQB,\r
  11   string &Q3, string &A3, string &B3);\r
  12 \r
  13 void AlignChimeLocal3(const string &Q3, const string &A3, const string &B3,\r
  14   const string &QLabel, const string &ALabel, const string &BLabel,\r
  15   ChimeHit2 &Hit);\r
  16 \r
  17 double GetScore2(double Y, double N, double A)\r
  18         {\r
  19         return Y/(opt_xn*(N + opt_dn) + opt_xa*A);\r
  20         }\r
  21 \r
  22 void AlignChimeGlobal3(const string &Q3, const string &A3, const string &B3,\r
  23   const string &QLabel, const string &ALabel, const string &BLabel,\r
  24   ChimeHit2 &Hit)\r
  25         {\r
  26         Hit.Clear();\r
  27         Hit.QLabel = QLabel;\r
  28 \r
  29         const byte *Q3Seq = (const byte *) Q3.c_str();\r
  30         const byte *A3Seq = (const byte *) A3.c_str();\r
  31         const byte *B3Seq = (const byte *) B3.c_str();\r
  32 \r
  33         const unsigned ColCount = SIZE(Q3);\r
  34         asserta(SIZE(A3) == ColCount && SIZE(B3) == ColCount);\r
  35 \r
  36 #if     TRACE\r
  37         Log("Q %5u %*.*s\n", ColCount, ColCount, ColCount, Q3Seq);\r
  38         Log("A %5u %*.*s\n", ColCount, ColCount, ColCount, A3Seq);\r
  39         Log("B %5u %*.*s\n", ColCount, ColCount, ColCount, B3Seq);\r
  40 #endif\r
  41 \r
  42 // Discard terminal gaps\r
  43         unsigned ColLo = UINT_MAX;\r
  44         unsigned ColHi = UINT_MAX;\r
  45         for (unsigned Col = 2; Col + 2 < ColCount; ++Col)\r
  46                 {\r
  47                 char q = Q3Seq[Col];\r
  48                 char a = A3Seq[Col];\r
  49                 char b = B3Seq[Col];\r
  50 \r
  51                 if (isacgt(q) && isacgt(a) && isacgt(b))\r
  52                         {\r
  53                         if (ColLo == UINT_MAX)\r
  54                                 ColLo = Col;\r
  55                         ColHi = Col;\r
  56                         }\r
  57                 }\r
  58 \r
  59         if (ColLo == UINT_MAX)\r
  60                 return;\r
  61 \r
  62         unsigned QPos = 0;\r
  63         unsigned APos = 0;\r
  64         unsigned BPos = 0;\r
  65         unsigned DiffCount = 0;\r
  66 \r
  67         vector<unsigned> ColToQPos(ColLo, UINT_MAX);\r
  68         vector<unsigned> AccumCount(ColLo, UINT_MAX);\r
  69         vector<unsigned> AccumSameA(ColLo, UINT_MAX);\r
  70         vector<unsigned> AccumSameB(ColLo, UINT_MAX);\r
  71         vector<unsigned> AccumForA(ColLo, UINT_MAX);\r
  72         vector<unsigned> AccumForB(ColLo, UINT_MAX);\r
  73         vector<unsigned> AccumAbstain(ColLo, UINT_MAX);\r
  74         vector<unsigned> AccumAgainst(ColLo, UINT_MAX);\r
  75 \r
  76         unsigned SumSameA = 0;\r
  77         unsigned SumSameB = 0;\r
  78         unsigned SumSameAB = 0;\r
  79         unsigned Sum = 0;\r
  80         unsigned SumForA = 0;\r
  81         unsigned SumForB = 0;\r
  82         unsigned SumAbstain = 0;\r
  83         unsigned SumAgainst = 0;\r
  84         for (unsigned Col = ColLo; Col <= ColHi; ++Col)\r
  85                 {\r
  86                 char q = Q3Seq[Col];\r
  87                 char a = A3Seq[Col];\r
  88                 char b = B3Seq[Col];\r
  89 \r
  90                 if (isacgt(q) && isacgt(a) && isacgt(b))\r
  91                         {\r
  92                         if (q == a)\r
  93                                 ++SumSameA;\r
  94                         if (q == b)\r
  95                                 ++SumSameB;\r
  96                         if (a == b)\r
  97                                 ++SumSameAB;\r
  98                         if (q == a && q != b)\r
  99                                 ++SumForA;\r
 100                         if (q == b && q != a)\r
 101                                 ++SumForB;\r
 102                         if (a == b && q != a)\r
 103                                 ++SumAgainst;\r
 104                         if (q != a && q != b)\r
 105                                 ++SumAbstain;\r
 106                         ++Sum;\r
 107                         }\r
 108 \r
 109                 ColToQPos.push_back(QPos);\r
 110                 AccumSameA.push_back(SumSameA);\r
 111                 AccumSameB.push_back(SumSameB);\r
 112                 AccumCount.push_back(Sum);\r
 113                 AccumForA.push_back(SumForA);\r
 114                 AccumForB.push_back(SumForB);\r
 115                 AccumAbstain.push_back(SumAbstain);\r
 116                 AccumAgainst.push_back(SumAgainst);\r
 117 \r
 118                 if (q != '-')\r
 119                         ++QPos;\r
 120                 if (a != '-')\r
 121                         ++APos;\r
 122                 if (b != '-')\r
 123                         ++BPos;\r
 124                 }\r
 125 \r
 126         asserta(SIZE(ColToQPos) == ColHi+1);\r
 127         asserta(SIZE(AccumSameA) == ColHi+1);\r
 128         asserta(SIZE(AccumSameB) == ColHi+1);\r
 129         asserta(SIZE(AccumAbstain) == ColHi+1);\r
 130         asserta(SIZE(AccumAgainst) == ColHi+1);\r
 131 \r
 132         double IdQA = double(SumSameA)/Sum;\r
 133         double IdQB = double(SumSameB)/Sum;\r
 134         double IdAB = double(SumSameAB)/Sum;\r
 135         double MaxId = max(IdQA, IdQB);\r
 136 \r
 137 #if     TRACE\r
 138         Log("IdQA=%.1f%% IdQB=%.1f%% IdAB=%.1f\n", IdQA*100.0, IdQB*100.0, IdAB*100.0);\r
 139         Log("\n");\r
 140         Log("    x  AQB   IdAL   IdBL   IdAR   IdBR   DivAB   DivBA    YAL    YBL    YAR    YBR    AbL    AbR  ScoreAB  ScoreAB    XLo    Xhi\n");\r
 141         Log("-----  ---  -----  -----  -----  -----  ------  ------  -----  -----  -----  -----  -----  -----  -------  -------  -----  -----\n");\r
 142 #endif\r
 143         unsigned BestXLo = UINT_MAX;\r
 144         unsigned BestXHi = UINT_MAX;\r
 145         double BestDiv = 0.0;\r
 146         double BestIdQM = 0.0;\r
 147         double BestScore = 0.0;\r
 148 \r
 149 // Find range of cols BestXLo..BestXHi that maximizes score\r
 150         bool FirstA = false;\r
 151 \r
 152 // NOTE: Must be < ColHi not <= because use Col+1 below\r
 153         for (unsigned Col = ColLo; Col < ColHi; ++Col)\r
 154                 {\r
 155                 char q = Q3Seq[Col];\r
 156                 char a = A3Seq[Col];\r
 157                 char b = B3Seq[Col];\r
 158 \r
 159                 unsigned SameAL = AccumSameA[Col];\r
 160                 unsigned SameBL = AccumSameB[Col];\r
 161                 unsigned SameAR = SumSameA - AccumSameA[Col];\r
 162                 unsigned SameBR = SumSameB - AccumSameB[Col];\r
 163 \r
 164                 double IdAB = double(SameAL + SameBR)/Sum;\r
 165                 double IdBA = double(SameBL + SameAR)/Sum;\r
 166 \r
 167                 unsigned ForAL = AccumForA[Col];\r
 168                 unsigned ForBL = AccumForB[Col];\r
 169                 unsigned ForAR = SumForA - AccumForA[Col+1];\r
 170                 unsigned ForBR = SumForB - AccumForB[Col+1];\r
 171                 unsigned AbL = AccumAbstain[Col];\r
 172                 unsigned AbR = SumAbstain - AccumAbstain[Col+1];\r
 173 \r
 174                 double ScoreAB = GetScore2(ForAL, ForBL, AbL)*GetScore2(ForBR, ForAR, AbR);\r
 175                 double ScoreBA = GetScore2(ForBL, ForAL, AbL)*GetScore2(ForAR, ForBR, AbR);\r
 176         \r
 177                 double DivAB = IdAB/MaxId;\r
 178                 double DivBA = IdBA/MaxId;\r
 179                 double MaxDiv = max(DivAB, DivBA);\r
 180 \r
 181                 //if (MaxDiv > BestDiv)\r
 182                 //      {\r
 183                 //      BestDiv = MaxDiv;\r
 184                 //      BestXLo = Col;\r
 185                 //      BestXHi = Col;\r
 186                 //      FirstA = (DivAB > DivBA);\r
 187                 //      if (FirstA)\r
 188                 //              BestIdQM = IdAB;\r
 189                 //      else\r
 190                 //              BestIdQM = IdBA;\r
 191                 //      }\r
 192                 //else if (MaxDiv == BestDiv)\r
 193                 //      BestXHi = Col;\r
 194 \r
 195                 double MaxScore = max(ScoreAB, ScoreBA);\r
 196                 if (MaxScore > BestScore)\r
 197                         {\r
 198                         BestScore = MaxScore;\r
 199                         BestXLo = Col;\r
 200                         BestXHi = Col;\r
 201                         FirstA = (ScoreAB > ScoreBA);\r
 202                         if (FirstA)\r
 203                                 BestIdQM = IdAB;\r
 204                         else\r
 205                                 BestIdQM = IdBA;\r
 206                         if (MaxDiv > BestDiv)\r
 207                                 BestDiv = MaxDiv;\r
 208                         }\r
 209                 else if (MaxScore == BestScore)\r
 210                         {\r
 211                         BestXHi = Col;\r
 212                         if (MaxDiv > BestDiv)\r
 213                                 BestDiv = MaxDiv;\r
 214                         }\r
 215 \r
 216 #if     TRACE\r
 217                 {\r
 218                 Log("%5u", Col);\r
 219                 char q = Q3Seq[Col];\r
 220                 char a = A3Seq[Col];\r
 221                 char b = B3Seq[Col];\r
 222                 Log("  %c%c%c", a, q, b);\r
 223                 Log("  %5u", SameAL);\r
 224                 Log("  %5u", SameBL);\r
 225                 Log("  %5u", SameAR);\r
 226                 Log("  %5u", SameBR);\r
 227                 Log("  %5.4f", DivAB);\r
 228                 Log("  %5.4f", DivBA);\r
 229                 Log("  %5u", ForAL);\r
 230                 Log("  %5u", ForBL);\r
 231                 Log("  %5u", ForAR);\r
 232                 Log("  %5u", ForBR);\r
 233                 Log("  %5u", AbL);\r
 234                 Log("  %5u", AbR);\r
 235                 Log("  %7.4f", ScoreAB);\r
 236                 Log("  %7.4f", ScoreBA);\r
 237                 if (BestXLo != UINT_MAX)\r
 238                         Log("  %5u", BestXLo);\r
 239                 if (BestXHi != UINT_MAX)\r
 240                         Log("  %5u", BestXHi);\r
 241                 Log("\n");\r
 242                 }\r
 243 #endif\r
 244                 }\r
 245 \r
 246         if (BestXLo == UINT_MAX)\r
 247                 {\r
 248 #if     TRACE\r
 249                 Log("\n");\r
 250                 Log("No crossover found.\n");\r
 251 #endif\r
 252                 return;\r
 253                 }\r
 254 #if     TRACE\r
 255         Log("BestX col %u - %u\n", BestXLo, BestXHi);\r
 256 #endif\r
 257 \r
 258 // Find maximum region of identity within BestXLo..BestXHi\r
 259         unsigned ColXLo = (BestXLo + BestXHi)/2;\r
 260         unsigned ColXHi = ColXLo;\r
 261         unsigned SegLo = UINT_MAX;\r
 262         unsigned SegHi = UINT_MAX;\r
 263         for (unsigned Col = BestXLo; Col <= BestXHi; ++Col)\r
 264                 {\r
 265                 char q = Q3Seq[Col];\r
 266                 char a = A3Seq[Col];\r
 267                 char b = B3Seq[Col];\r
 268 \r
 269                 if (q == a && q == b)\r
 270                         {\r
 271                         if (SegLo == UINT_MAX)\r
 272                                 SegLo = Col;\r
 273                         SegHi = Col;\r
 274                         }\r
 275                 else\r
 276                         {\r
 277                         unsigned SegLength = SegHi - SegLo + 1;\r
 278                         unsigned BestSegLength = ColXHi - ColXLo + 1;\r
 279                         if (SegLength > BestSegLength)\r
 280                                 {\r
 281                                 ColXLo = SegLo;\r
 282                                 ColXHi = SegHi;\r
 283                                 }\r
 284                         SegLo = UINT_MAX;\r
 285                         SegHi = UINT_MAX;\r
 286                         }\r
 287                 }\r
 288         unsigned SegLength = SegHi - SegLo + 1;\r
 289         unsigned BestSegLength = ColXHi - ColXLo + 1;\r
 290         if (SegLength > BestSegLength)\r
 291                 {\r
 292                 ColXLo = SegLo;\r
 293                 ColXHi = SegHi;\r
 294                 }\r
 295 \r
 296         QPos = 0;\r
 297         for (unsigned x = 0; x < ColCount; ++x)\r
 298                 {\r
 299                 if (x == ColXLo)\r
 300                         Hit.QXLo = QPos;\r
 301                 else if (x == ColXHi)\r
 302                         {\r
 303                         Hit.QXHi = QPos;\r
 304                         break;\r
 305                         }\r
 306                 char q = Q3Seq[x];\r
 307                 if (q != '-')\r
 308                         ++QPos;\r
 309                 }\r
 310 \r
 311         Hit.ColXLo = ColXLo;\r
 312         Hit.ColXHi = ColXHi;\r
 313 \r
 314         //if (FirstA)\r
 315         //      {\r
 316         //      Hit.LY = AccumForA[ColXLo];\r
 317         //      Hit.LN = AccumForB[ColXLo];\r
 318 \r
 319         //      Hit.RY = SumForB - AccumForB[ColXHi];\r
 320         //      Hit.RN = SumForA - AccumForA[ColXHi];\r
 321         //      }\r
 322         //else\r
 323         //      {\r
 324         //      Hit.LY = AccumForB[ColXLo];\r
 325         //      Hit.LN = AccumForA[ColXLo];\r
 326         //      Hit.RY = SumForA - AccumForA[ColXHi];\r
 327         //      Hit.RN = SumForB - AccumForB[ColXHi];\r
 328         //      }\r
 329 \r
 330         //Hit.LA = AccumAgainst[ColXLo];\r
 331         //Hit.LD = AccumAbstain[ColXLo];\r
 332 \r
 333         //Hit.RA = SumAgainst - AccumAgainst[ColXHi];\r
 334         //Hit.RD = SumAbstain - AccumAbstain[ColXHi];\r
 335 \r
 336         Hit.PctIdAB = IdAB*100.0;\r
 337         Hit.PctIdQM = BestIdQM*100.0;\r
 338 \r
 339         Hit.Div = (BestDiv - 1.0)*100.0;\r
 340 \r
 341         //Hit.QSD = QSD;\r
 342         Hit.Q3 = Q3;\r
 343         Hit.QLabel = QLabel;\r
 344         if (FirstA)\r
 345                 {\r
 346                 //Hit.ASD = ASD;\r
 347                 //Hit.BSD = BSD;\r
 348                 //Hit.PathQA = PathQA;\r
 349                 //Hit.PathQB = PathQB;\r
 350                 Hit.A3 = A3;\r
 351                 Hit.B3 = B3;\r
 352                 Hit.ALabel = ALabel;\r
 353                 Hit.BLabel = BLabel;\r
 354                 Hit.PctIdQA = IdQA*100.0;\r
 355                 Hit.PctIdQB = IdQB*100.0;\r
 356                 }\r
 357         else\r
 358                 {\r
 359                 Hit.A3 = B3;\r
 360                 Hit.B3 = A3;\r
 361                 Hit.ALabel = BLabel;\r
 362                 Hit.BLabel = ALabel;\r
 363                 Hit.PctIdQA = IdQB*100.0;\r
 364                 Hit.PctIdQB = IdQA*100.0;\r
 365                 }\r
 366 \r
 367 // CS SNPs\r
 368         Hit.CS_LY = 0;\r
 369         Hit.CS_LN = 0;\r
 370         Hit.CS_RY = 0;\r
 371         Hit.CS_RN = 0;\r
 372         Hit.CS_LA = 0;\r
 373         Hit.CS_RA = 0;\r
 374 \r
 375         //vector<float> Cons;\r
 376         //for (unsigned Col = 0; Col < ColCount; ++Col)\r
 377         //      {\r
 378         //      char q = Q3Seq[Col];\r
 379         //      char a = A3Seq[Col];\r
 380         //      char b = B3Seq[Col];\r
 381         //      if (q == a && q == b && a == b)\r
 382         //              {\r
 383         //              Cons.push_back(1.0f);\r
 384         //              continue;\r
 385         //              }\r
 386 \r
 387         //      bool gapq = isgap(q);\r
 388         //      bool gapa = isgap(a);\r
 389         //      bool gapb = isgap(b);\r
 390 \r
 391         //      if (!gapq && !gapa && !gapb)\r
 392         //              {\r
 393         //              if (q == a || q == b || a == b)\r
 394         //                      Cons.push_back(0.75);\r
 395         //              else\r
 396         //                      Cons.push_back(0.5);\r
 397         //              }\r
 398         //      else\r
 399         //              {\r
 400         //              if (!gapa && (a == b || a == q))\r
 401         //                      Cons.push_back(0.5f);\r
 402         //              else if (!gapb && b == q)\r
 403         //                      Cons.push_back(0.5f);\r
 404         //              else\r
 405         //                      Cons.push_back(0.0f);\r
 406         //              }\r
 407         //      }\r
 408 \r
 409         //float fLY = 0.0f;\r
 410         //float fLN = 0.0f;\r
 411         //float fLA = 0.0f;\r
 412         //float fRY = 0.0f;\r
 413         //float fRN = 0.0f;\r
 414         //float fRA = 0.0f;\r
 415         for (unsigned Col = ColLo; Col <= ColHi; ++Col)\r
 416                 {\r
 417                 char q = Q3Seq[Col];\r
 418                 char a = A3Seq[Col];\r
 419                 char b = B3Seq[Col];\r
 420                 if (q == a && q == b && a == b)\r
 421                         continue;\r
 422 \r
 423                 unsigned ngaps = 0;\r
 424                 if (isgap(q))\r
 425                         ++ngaps;\r
 426                 if (isgap(a))\r
 427                         ++ngaps;\r
 428                 if (isgap(b))\r
 429                         ++ngaps;\r
 430 \r
 431                 if (opt_skipgaps)\r
 432                         {\r
 433                         if (ngaps == 3)\r
 434                                 continue;\r
 435                         }\r
 436                 else\r
 437                         {\r
 438                         if (ngaps == 2)\r
 439                                 continue;\r
 440                         }\r
 441 \r
 442                 if (!FirstA)\r
 443                         swap(a, b);\r
 444 \r
 445                 //float AvgCons = (Cons[Col-2] + Cons[Col-1] + Cons[Col+1] + Cons[Col+2])/4;\r
 446                 //if (Col < ColXLo)\r
 447                 //      {\r
 448                 //      if (q == a && q != b)\r
 449                 //              fLY += AvgCons;\r
 450                 //      else if (q == b && q != a)\r
 451                 //              fLN += AvgCons;\r
 452                 //      else\r
 453                 //              fLA += AvgCons;\r
 454                 //      }\r
 455                 //else if (Col > ColXHi)\r
 456                 //      {\r
 457                 //      if (q == b && q != a)\r
 458                 //              fRY += AvgCons;\r
 459                 //      else if (q == a && q != b)\r
 460                 //              fRN += AvgCons;\r
 461                 //      else\r
 462                 //              fRA += AvgCons;\r
 463                 //      }\r
 464 \r
 465                 if (opt_skipgaps2)\r
 466                         {\r
 467                         if (Col > 0 && (isgap(Q3Seq[Col-1]) || isgap(A3Seq[Col-1]) || isgap(B3Seq[Col-1])))\r
 468                                 continue;\r
 469                         if (Col + 1 < ColCount && (isgap(Q3Seq[Col+1]) || isgap(A3Seq[Col+1]) || isgap(B3Seq[Col+1])))\r
 470                                 continue;\r
 471                         }\r
 472 \r
 473                 //if (Col > 0 && isgap(Q3Seq[Col-1]))\r
 474                         //continue;\r
 475                 //if (Col + 1 < ColCount && isgap(Q3Seq[Col+1]))\r
 476                 //      continue;\r
 477 \r
 478                 if (Col < ColXLo)\r
 479                         {\r
 480                         if (q == a && q != b)\r
 481                                 ++Hit.CS_LY;\r
 482                         else if (q == b && q != a)\r
 483                                 ++Hit.CS_LN;\r
 484                         else\r
 485                                 ++Hit.CS_LA;\r
 486                         }\r
 487                 else if (Col > ColXHi)\r
 488                         {\r
 489                         if (q == b && q != a)\r
 490                                 ++Hit.CS_RY;\r
 491                         else if (q == a && q != b)\r
 492                                 ++Hit.CS_RN;\r
 493                         else\r
 494                                 ++Hit.CS_RA;\r
 495                         }\r
 496                 }\r
 497 \r
 498         double ScoreL = GetScore2(Hit.CS_LY, Hit.CS_LN, Hit.CS_LA);\r
 499         double ScoreR = GetScore2(Hit.CS_RY, Hit.CS_RN, Hit.CS_RA);\r
 500         Hit.Score = ScoreL*ScoreR;\r
 501 \r
 502         extern bool g_UchimeDeNovo;\r
 503 \r
 504         //if (0)//g_UchimeDeNovo)\r
 505         //      {\r
 506         //      double AbQ = GetAbFromLabel(QLabel.c_str());\r
 507         //      double AbA = GetAbFromLabel(ALabel.c_str());\r
 508         //      double AbB = GetAbFromLabel(BLabel.c_str());\r
 509         //      if (AbQ > 0.0 && AbA > 0.0 && AbB > 0.0)\r
 510         //              {\r
 511         //              double MinAb = min(AbA, AbB);\r
 512         //              double Ratio = MinAb/AbQ;\r
 513         //              double t = Ratio - opt_abx;\r
 514         //      //      double Factor = 2.0/(1.0 + exp(-t));\r
 515         //              double Factor = min(Ratio, opt_abx)/opt_abx;\r
 516         //              if (opt_verbose)\r
 517         //                      Log("Score %.4f Ab factor %.4f >%s\n", Hit.Score, Factor, QLabel.c_str());\r
 518         //              Hit.Score *= Factor;\r
 519         //              }\r
 520         //      }\r
 521 \r
 522         extern FILE *g_fUChimeAlns;\r
 523         if (g_fUChimeAlns != 0 && Hit.Div > 0.0)\r
 524                 {\r
 525                 void WriteChimeHitX(FILE *f, const ChimeHit2 &Hit);\r
 526                 WriteChimeHitX(g_fUChimeAlns, Hit);\r
 527                 }\r
 528         }\r
 529 \r
 530 void AlignChime3(const string &Q3, const string &A3, const string &B3,\r
 531   const string &QLabel, const string &ALabel, const string &BLabel,\r
 532   ChimeHit2 &Hit)\r
 533         {\r
 534         if (opt_ucl)\r
 535                 AlignChimeLocal3(Q3, A3, B3, QLabel, ALabel, BLabel, Hit);\r
 536         else\r
 537                 AlignChimeGlobal3(Q3, A3, B3, QLabel, ALabel, BLabel, Hit);\r
 538         }\r
 539 \r
 540 static void StripGaps(const byte *Seq, unsigned L, string &s)\r
 541         {\r
 542         s.clear();\r
 543         for (unsigned i = 0; i < L; ++i)\r
 544                 {\r
 545                 char c = Seq[i];\r
 546                 if (!isgap(c))\r
 547                         s.push_back(c);\r
 548                 }\r
 549         }\r
 550 \r
 551 static void StripGapsAlloc(const SeqData &SDIn, SeqData &SDOut)\r
 552         {\r
 553         SDOut = SDIn;\r
 554         byte *s = myalloc(byte, SDIn.L);\r
 555         unsigned k = 0;\r
 556         for (unsigned i = 0; i < SDIn.L; ++i)\r
 557                 {\r
 558                 char c = SDIn.Seq[i];\r
 559                 if (!isgap(c))\r
 560                         s[k++] = toupper(c);\r
 561                 }\r
 562         SDOut.Seq = s;\r
 563         SDOut.L = k;\r
 564         }\r
 565 \r
 566 void AlignChime(const SeqData &QSD, const SeqData &ASD, const SeqData &BSD,\r
 567   const string &PathQA, const string &PathQB, ChimeHit2 &Hit)\r
 568         {\r
 569         //if (opt_ucl)\r
 570         //      {\r
 571         //      AlignChimeLocal(QSD, ASD, BSD, PathQA, PathQB, Hit);\r
 572         //      return;\r
 573         //      }\r
 574 \r
 575         string Q3;\r
 576         string A3;\r
 577         string B3;\r
 578         Make3Way(QSD, ASD, BSD, PathQA, PathQB, Q3, A3, B3);\r
 579 \r
 580         AlignChime3(Q3, A3, B3, QSD.Label, ASD.Label, BSD.Label, Hit);\r
 581         }\r
 582 \r
 583 void AlignChime3SDRealign(const SeqData &QSD3, const SeqData &ASD3, const SeqData &BSD3,\r
 584   ChimeHit2 &Hit)\r
 585         {\r
 586         SeqData QSD;\r
 587         SeqData ASD;\r
 588         SeqData BSD;\r
 589         StripGapsAlloc(QSD3, QSD);\r
 590         StripGapsAlloc(ASD3, ASD);\r
 591         StripGapsAlloc(BSD3, BSD);\r
 592 \r
 593         string PathQA;\r
 594         string PathQB;\r
 595         bool FoundQA = GlobalAlign(QSD, ASD, PathQA);\r
 596         bool FoundQB = GlobalAlign(QSD, BSD, PathQB);\r
 597         if (!FoundQA || !FoundQB)\r
 598                 {\r
 599                 Hit.Clear();\r
 600                 Hit.QLabel = QSD3.Label;\r
 601                 return;\r
 602                 }\r
 603 \r
 604         AlignChime(QSD, ASD, BSD, PathQA, PathQB, Hit);\r
 605 \r
 606         myfree((void *) QSD.Seq);\r
 607         myfree((void *) ASD.Seq);\r
 608         myfree((void *) BSD.Seq);\r
 609         }\r
 610 \r
 611 void AlignChime3SD(const SeqData &QSD3, const SeqData &ASD3, const SeqData &BSD3,\r
 612   ChimeHit2 &Hit)\r
 613         {\r
 614         if (opt_realign)\r
 615                 {\r
 616                 AlignChime3SDRealign(QSD3, ASD3, BSD3, Hit);\r
 617                 return;\r
 618                 }\r
 619 \r
 620         string Q3;\r
 621         string A3;\r
 622         string B3;\r
 623 \r
 624         const unsigned ColCount = QSD3.L;\r
 625         asserta(ASD3.L == ColCount && BSD3.L == ColCount);\r
 626 \r
 627         Q3.reserve(ColCount);\r
 628         A3.reserve(ColCount);\r
 629         B3.reserve(ColCount);\r
 630 \r
 631         const byte *QS = QSD3.Seq;\r
 632         const byte *AS = ASD3.Seq;\r
 633         const byte *BS = BSD3.Seq;\r
 634         for (unsigned Col = 0; Col < ColCount; ++Col)\r
 635                 {\r
 636                 byte q = toupper(QS[Col]);\r
 637                 byte a = toupper(AS[Col]);\r
 638                 byte b = toupper(BS[Col]);\r
 639 \r
 640                 if (isgap(q) && isgap(a) && isgap(b))\r
 641                         continue;\r
 642 \r
 643                 Q3.push_back(q);\r
 644                 A3.push_back(a);\r
 645                 B3.push_back(b);\r
 646                 }\r
 647 \r
 648         AlignChime3(Q3, A3, B3, QSD3.Label, ASD3.Label, BSD3.Label, Hit);\r
 649         }\r