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