--- /dev/null
+#include "myutils.h"\r
+#include "seq.h"\r
+#include "chime.h"\r
+#include "dp.h"\r
+\r
+#define TRACE 0\r
+#define TRACE_BS 0\r
+\r
+void Make3Way(const SeqData &SDQ, const SeqData &SDA, const SeqData &SDB,\r
+ const string &PathQA, const string &PathQB,\r
+ string &Q3, string &A3, string &B3);\r
+\r
+void AlignChimeLocal3(const string &Q3, const string &A3, const string &B3,\r
+ const string &QLabel, const string &ALabel, const string &BLabel,\r
+ ChimeHit2 &Hit);\r
+\r
+double GetScore2(double Y, double N, double A)\r
+ {\r
+ return Y/(opt_xn*(N + opt_dn) + opt_xa*A);\r
+ }\r
+\r
+void AlignChimeGlobal3(const string &Q3, const string &A3, const string &B3,\r
+ const string &QLabel, const string &ALabel, const string &BLabel,\r
+ ChimeHit2 &Hit)\r
+ {\r
+ Hit.Clear();\r
+ Hit.QLabel = QLabel;\r
+\r
+ const byte *Q3Seq = (const byte *) Q3.c_str();\r
+ const byte *A3Seq = (const byte *) A3.c_str();\r
+ const byte *B3Seq = (const byte *) B3.c_str();\r
+\r
+ const unsigned ColCount = SIZE(Q3);\r
+ asserta(SIZE(A3) == ColCount && SIZE(B3) == ColCount);\r
+\r
+#if TRACE\r
+ Log("Q %5u %*.*s\n", ColCount, ColCount, ColCount, Q3Seq);\r
+ Log("A %5u %*.*s\n", ColCount, ColCount, ColCount, A3Seq);\r
+ Log("B %5u %*.*s\n", ColCount, ColCount, ColCount, B3Seq);\r
+#endif\r
+\r
+// Discard terminal gaps\r
+ unsigned ColLo = UINT_MAX;\r
+ unsigned ColHi = UINT_MAX;\r
+ for (unsigned Col = 2; Col + 2 < ColCount; ++Col)\r
+ {\r
+ char q = Q3Seq[Col];\r
+ char a = A3Seq[Col];\r
+ char b = B3Seq[Col];\r
+\r
+ if (isacgt(q) && isacgt(a) && isacgt(b))\r
+ {\r
+ if (ColLo == UINT_MAX)\r
+ ColLo = Col;\r
+ ColHi = Col;\r
+ }\r
+ }\r
+\r
+ if (ColLo == UINT_MAX)\r
+ return;\r
+\r
+ unsigned QPos = 0;\r
+ unsigned APos = 0;\r
+ unsigned BPos = 0;\r
+ unsigned DiffCount = 0;\r
+\r
+ vector<unsigned> ColToQPos(ColLo, UINT_MAX);\r
+ vector<unsigned> AccumCount(ColLo, UINT_MAX);\r
+ vector<unsigned> AccumSameA(ColLo, UINT_MAX);\r
+ vector<unsigned> AccumSameB(ColLo, UINT_MAX);\r
+ vector<unsigned> AccumForA(ColLo, UINT_MAX);\r
+ vector<unsigned> AccumForB(ColLo, UINT_MAX);\r
+ vector<unsigned> AccumAbstain(ColLo, UINT_MAX);\r
+ vector<unsigned> AccumAgainst(ColLo, UINT_MAX);\r
+\r
+ unsigned SumSameA = 0;\r
+ unsigned SumSameB = 0;\r
+ unsigned SumSameAB = 0;\r
+ unsigned Sum = 0;\r
+ unsigned SumForA = 0;\r
+ unsigned SumForB = 0;\r
+ unsigned SumAbstain = 0;\r
+ unsigned SumAgainst = 0;\r
+ for (unsigned Col = ColLo; Col <= ColHi; ++Col)\r
+ {\r
+ char q = Q3Seq[Col];\r
+ char a = A3Seq[Col];\r
+ char b = B3Seq[Col];\r
+\r
+ if (isacgt(q) && isacgt(a) && isacgt(b))\r
+ {\r
+ if (q == a)\r
+ ++SumSameA;\r
+ if (q == b)\r
+ ++SumSameB;\r
+ if (a == b)\r
+ ++SumSameAB;\r
+ if (q == a && q != b)\r
+ ++SumForA;\r
+ if (q == b && q != a)\r
+ ++SumForB;\r
+ if (a == b && q != a)\r
+ ++SumAgainst;\r
+ if (q != a && q != b)\r
+ ++SumAbstain;\r
+ ++Sum;\r
+ }\r
+\r
+ ColToQPos.push_back(QPos);\r
+ AccumSameA.push_back(SumSameA);\r
+ AccumSameB.push_back(SumSameB);\r
+ AccumCount.push_back(Sum);\r
+ AccumForA.push_back(SumForA);\r
+ AccumForB.push_back(SumForB);\r
+ AccumAbstain.push_back(SumAbstain);\r
+ AccumAgainst.push_back(SumAgainst);\r
+\r
+ if (q != '-')\r
+ ++QPos;\r
+ if (a != '-')\r
+ ++APos;\r
+ if (b != '-')\r
+ ++BPos;\r
+ }\r
+\r
+ asserta(SIZE(ColToQPos) == ColHi+1);\r
+ asserta(SIZE(AccumSameA) == ColHi+1);\r
+ asserta(SIZE(AccumSameB) == ColHi+1);\r
+ asserta(SIZE(AccumAbstain) == ColHi+1);\r
+ asserta(SIZE(AccumAgainst) == ColHi+1);\r
+\r
+ double IdQA = double(SumSameA)/Sum;\r
+ double IdQB = double(SumSameB)/Sum;\r
+ double IdAB = double(SumSameAB)/Sum;\r
+ double MaxId = max(IdQA, IdQB);\r
+\r
+#if TRACE\r
+ Log("IdQA=%.1f%% IdQB=%.1f%% IdAB=%.1f\n", IdQA*100.0, IdQB*100.0, IdAB*100.0);\r
+ Log("\n");\r
+ Log(" x AQB IdAL IdBL IdAR IdBR DivAB DivBA YAL YBL YAR YBR AbL AbR ScoreAB ScoreAB XLo Xhi\n");\r
+ Log("----- --- ----- ----- ----- ----- ------ ------ ----- ----- ----- ----- ----- ----- ------- ------- ----- -----\n");\r
+#endif\r
+ unsigned BestXLo = UINT_MAX;\r
+ unsigned BestXHi = UINT_MAX;\r
+ double BestDiv = 0.0;\r
+ double BestIdQM = 0.0;\r
+ double BestScore = 0.0;\r
+\r
+// Find range of cols BestXLo..BestXHi that maximizes score\r
+ bool FirstA = false;\r
+\r
+// NOTE: Must be < ColHi not <= because use Col+1 below\r
+ for (unsigned Col = ColLo; Col < ColHi; ++Col)\r
+ {\r
+ char q = Q3Seq[Col];\r
+ char a = A3Seq[Col];\r
+ char b = B3Seq[Col];\r
+\r
+ unsigned SameAL = AccumSameA[Col];\r
+ unsigned SameBL = AccumSameB[Col];\r
+ unsigned SameAR = SumSameA - AccumSameA[Col];\r
+ unsigned SameBR = SumSameB - AccumSameB[Col];\r
+\r
+ double IdAB = double(SameAL + SameBR)/Sum;\r
+ double IdBA = double(SameBL + SameAR)/Sum;\r
+\r
+ unsigned ForAL = AccumForA[Col];\r
+ unsigned ForBL = AccumForB[Col];\r
+ unsigned ForAR = SumForA - AccumForA[Col+1];\r
+ unsigned ForBR = SumForB - AccumForB[Col+1];\r
+ unsigned AbL = AccumAbstain[Col];\r
+ unsigned AbR = SumAbstain - AccumAbstain[Col+1];\r
+\r
+ double ScoreAB = GetScore2(ForAL, ForBL, AbL)*GetScore2(ForBR, ForAR, AbR);\r
+ double ScoreBA = GetScore2(ForBL, ForAL, AbL)*GetScore2(ForAR, ForBR, AbR);\r
+ \r
+ double DivAB = IdAB/MaxId;\r
+ double DivBA = IdBA/MaxId;\r
+ double MaxDiv = max(DivAB, DivBA);\r
+\r
+ //if (MaxDiv > BestDiv)\r
+ // {\r
+ // BestDiv = MaxDiv;\r
+ // BestXLo = Col;\r
+ // BestXHi = Col;\r
+ // FirstA = (DivAB > DivBA);\r
+ // if (FirstA)\r
+ // BestIdQM = IdAB;\r
+ // else\r
+ // BestIdQM = IdBA;\r
+ // }\r
+ //else if (MaxDiv == BestDiv)\r
+ // BestXHi = Col;\r
+\r
+ double MaxScore = max(ScoreAB, ScoreBA);\r
+ if (MaxScore > BestScore)\r
+ {\r
+ BestScore = MaxScore;\r
+ BestXLo = Col;\r
+ BestXHi = Col;\r
+ FirstA = (ScoreAB > ScoreBA);\r
+ if (FirstA)\r
+ BestIdQM = IdAB;\r
+ else\r
+ BestIdQM = IdBA;\r
+ if (MaxDiv > BestDiv)\r
+ BestDiv = MaxDiv;\r
+ }\r
+ else if (MaxScore == BestScore)\r
+ {\r
+ BestXHi = Col;\r
+ if (MaxDiv > BestDiv)\r
+ BestDiv = MaxDiv;\r
+ }\r
+\r
+#if TRACE\r
+ {\r
+ Log("%5u", Col);\r
+ char q = Q3Seq[Col];\r
+ char a = A3Seq[Col];\r
+ char b = B3Seq[Col];\r
+ Log(" %c%c%c", a, q, b);\r
+ Log(" %5u", SameAL);\r
+ Log(" %5u", SameBL);\r
+ Log(" %5u", SameAR);\r
+ Log(" %5u", SameBR);\r
+ Log(" %5.4f", DivAB);\r
+ Log(" %5.4f", DivBA);\r
+ Log(" %5u", ForAL);\r
+ Log(" %5u", ForBL);\r
+ Log(" %5u", ForAR);\r
+ Log(" %5u", ForBR);\r
+ Log(" %5u", AbL);\r
+ Log(" %5u", AbR);\r
+ Log(" %7.4f", ScoreAB);\r
+ Log(" %7.4f", ScoreBA);\r
+ if (BestXLo != UINT_MAX)\r
+ Log(" %5u", BestXLo);\r
+ if (BestXHi != UINT_MAX)\r
+ Log(" %5u", BestXHi);\r
+ Log("\n");\r
+ }\r
+#endif\r
+ }\r
+\r
+ if (BestXLo == UINT_MAX)\r
+ {\r
+#if TRACE\r
+ Log("\n");\r
+ Log("No crossover found.\n");\r
+#endif\r
+ return;\r
+ }\r
+#if TRACE\r
+ Log("BestX col %u - %u\n", BestXLo, BestXHi);\r
+#endif\r
+\r
+// Find maximum region of identity within BestXLo..BestXHi\r
+ unsigned ColXLo = (BestXLo + BestXHi)/2;\r
+ unsigned ColXHi = ColXLo;\r
+ unsigned SegLo = UINT_MAX;\r
+ unsigned SegHi = UINT_MAX;\r
+ for (unsigned Col = BestXLo; Col <= BestXHi; ++Col)\r
+ {\r
+ char q = Q3Seq[Col];\r
+ char a = A3Seq[Col];\r
+ char b = B3Seq[Col];\r
+\r
+ if (q == a && q == b)\r
+ {\r
+ if (SegLo == UINT_MAX)\r
+ SegLo = Col;\r
+ SegHi = Col;\r
+ }\r
+ else\r
+ {\r
+ unsigned SegLength = SegHi - SegLo + 1;\r
+ unsigned BestSegLength = ColXHi - ColXLo + 1;\r
+ if (SegLength > BestSegLength)\r
+ {\r
+ ColXLo = SegLo;\r
+ ColXHi = SegHi;\r
+ }\r
+ SegLo = UINT_MAX;\r
+ SegHi = UINT_MAX;\r
+ }\r
+ }\r
+ unsigned SegLength = SegHi - SegLo + 1;\r
+ unsigned BestSegLength = ColXHi - ColXLo + 1;\r
+ if (SegLength > BestSegLength)\r
+ {\r
+ ColXLo = SegLo;\r
+ ColXHi = SegHi;\r
+ }\r
+\r
+ QPos = 0;\r
+ for (unsigned x = 0; x < ColCount; ++x)\r
+ {\r
+ if (x == ColXLo)\r
+ Hit.QXLo = QPos;\r
+ else if (x == ColXHi)\r
+ {\r
+ Hit.QXHi = QPos;\r
+ break;\r
+ }\r
+ char q = Q3Seq[x];\r
+ if (q != '-')\r
+ ++QPos;\r
+ }\r
+\r
+ Hit.ColXLo = ColXLo;\r
+ Hit.ColXHi = ColXHi;\r
+\r
+ //if (FirstA)\r
+ // {\r
+ // Hit.LY = AccumForA[ColXLo];\r
+ // Hit.LN = AccumForB[ColXLo];\r
+\r
+ // Hit.RY = SumForB - AccumForB[ColXHi];\r
+ // Hit.RN = SumForA - AccumForA[ColXHi];\r
+ // }\r
+ //else\r
+ // {\r
+ // Hit.LY = AccumForB[ColXLo];\r
+ // Hit.LN = AccumForA[ColXLo];\r
+ // Hit.RY = SumForA - AccumForA[ColXHi];\r
+ // Hit.RN = SumForB - AccumForB[ColXHi];\r
+ // }\r
+\r
+ //Hit.LA = AccumAgainst[ColXLo];\r
+ //Hit.LD = AccumAbstain[ColXLo];\r
+\r
+ //Hit.RA = SumAgainst - AccumAgainst[ColXHi];\r
+ //Hit.RD = SumAbstain - AccumAbstain[ColXHi];\r
+\r
+ Hit.PctIdAB = IdAB*100.0;\r
+ Hit.PctIdQM = BestIdQM*100.0;\r
+\r
+ Hit.Div = (BestDiv - 1.0)*100.0;\r
+\r
+ //Hit.QSD = QSD;\r
+ Hit.Q3 = Q3;\r
+ Hit.QLabel = QLabel;\r
+ if (FirstA)\r
+ {\r
+ //Hit.ASD = ASD;\r
+ //Hit.BSD = BSD;\r
+ //Hit.PathQA = PathQA;\r
+ //Hit.PathQB = PathQB;\r
+ Hit.A3 = A3;\r
+ Hit.B3 = B3;\r
+ Hit.ALabel = ALabel;\r
+ Hit.BLabel = BLabel;\r
+ Hit.PctIdQA = IdQA*100.0;\r
+ Hit.PctIdQB = IdQB*100.0;\r
+ }\r
+ else\r
+ {\r
+ Hit.A3 = B3;\r
+ Hit.B3 = A3;\r
+ Hit.ALabel = BLabel;\r
+ Hit.BLabel = ALabel;\r
+ Hit.PctIdQA = IdQB*100.0;\r
+ Hit.PctIdQB = IdQA*100.0;\r
+ }\r
+\r
+// CS SNPs\r
+ Hit.CS_LY = 0;\r
+ Hit.CS_LN = 0;\r
+ Hit.CS_RY = 0;\r
+ Hit.CS_RN = 0;\r
+ Hit.CS_LA = 0;\r
+ Hit.CS_RA = 0;\r
+\r
+ //vector<float> Cons;\r
+ //for (unsigned Col = 0; Col < ColCount; ++Col)\r
+ // {\r
+ // char q = Q3Seq[Col];\r
+ // char a = A3Seq[Col];\r
+ // char b = B3Seq[Col];\r
+ // if (q == a && q == b && a == b)\r
+ // {\r
+ // Cons.push_back(1.0f);\r
+ // continue;\r
+ // }\r
+\r
+ // bool gapq = isgap(q);\r
+ // bool gapa = isgap(a);\r
+ // bool gapb = isgap(b);\r
+\r
+ // if (!gapq && !gapa && !gapb)\r
+ // {\r
+ // if (q == a || q == b || a == b)\r
+ // Cons.push_back(0.75);\r
+ // else\r
+ // Cons.push_back(0.5);\r
+ // }\r
+ // else\r
+ // {\r
+ // if (!gapa && (a == b || a == q))\r
+ // Cons.push_back(0.5f);\r
+ // else if (!gapb && b == q)\r
+ // Cons.push_back(0.5f);\r
+ // else\r
+ // Cons.push_back(0.0f);\r
+ // }\r
+ // }\r
+\r
+ //float fLY = 0.0f;\r
+ //float fLN = 0.0f;\r
+ //float fLA = 0.0f;\r
+ //float fRY = 0.0f;\r
+ //float fRN = 0.0f;\r
+ //float fRA = 0.0f;\r
+ for (unsigned Col = ColLo; Col <= ColHi; ++Col)\r
+ {\r
+ char q = Q3Seq[Col];\r
+ char a = A3Seq[Col];\r
+ char b = B3Seq[Col];\r
+ if (q == a && q == b && a == b)\r
+ continue;\r
+\r
+ unsigned ngaps = 0;\r
+ if (isgap(q))\r
+ ++ngaps;\r
+ if (isgap(a))\r
+ ++ngaps;\r
+ if (isgap(b))\r
+ ++ngaps;\r
+\r
+ if (opt_skipgaps)\r
+ {\r
+ if (ngaps == 3)\r
+ continue;\r
+ }\r
+ else\r
+ {\r
+ if (ngaps == 2)\r
+ continue;\r
+ }\r
+\r
+ if (!FirstA)\r
+ swap(a, b);\r
+\r
+ //float AvgCons = (Cons[Col-2] + Cons[Col-1] + Cons[Col+1] + Cons[Col+2])/4;\r
+ //if (Col < ColXLo)\r
+ // {\r
+ // if (q == a && q != b)\r
+ // fLY += AvgCons;\r
+ // else if (q == b && q != a)\r
+ // fLN += AvgCons;\r
+ // else\r
+ // fLA += AvgCons;\r
+ // }\r
+ //else if (Col > ColXHi)\r
+ // {\r
+ // if (q == b && q != a)\r
+ // fRY += AvgCons;\r
+ // else if (q == a && q != b)\r
+ // fRN += AvgCons;\r
+ // else\r
+ // fRA += AvgCons;\r
+ // }\r
+\r
+ if (opt_skipgaps2)\r
+ {\r
+ if (Col > 0 && (isgap(Q3Seq[Col-1]) || isgap(A3Seq[Col-1]) || isgap(B3Seq[Col-1])))\r
+ continue;\r
+ if (Col + 1 < ColCount && (isgap(Q3Seq[Col+1]) || isgap(A3Seq[Col+1]) || isgap(B3Seq[Col+1])))\r
+ continue;\r
+ }\r
+\r
+ //if (Col > 0 && isgap(Q3Seq[Col-1]))\r
+ //continue;\r
+ //if (Col + 1 < ColCount && isgap(Q3Seq[Col+1]))\r
+ // continue;\r
+\r
+ if (Col < ColXLo)\r
+ {\r
+ if (q == a && q != b)\r
+ ++Hit.CS_LY;\r
+ else if (q == b && q != a)\r
+ ++Hit.CS_LN;\r
+ else\r
+ ++Hit.CS_LA;\r
+ }\r
+ else if (Col > ColXHi)\r
+ {\r
+ if (q == b && q != a)\r
+ ++Hit.CS_RY;\r
+ else if (q == a && q != b)\r
+ ++Hit.CS_RN;\r
+ else\r
+ ++Hit.CS_RA;\r
+ }\r
+ }\r
+\r
+ double ScoreL = GetScore2(Hit.CS_LY, Hit.CS_LN, Hit.CS_LA);\r
+ double ScoreR = GetScore2(Hit.CS_RY, Hit.CS_RN, Hit.CS_RA);\r
+ Hit.Score = ScoreL*ScoreR;\r
+\r
+ extern bool g_UchimeDeNovo;\r
+\r
+ //if (0)//g_UchimeDeNovo)\r
+ // {\r
+ // double AbQ = GetAbFromLabel(QLabel.c_str());\r
+ // double AbA = GetAbFromLabel(ALabel.c_str());\r
+ // double AbB = GetAbFromLabel(BLabel.c_str());\r
+ // if (AbQ > 0.0 && AbA > 0.0 && AbB > 0.0)\r
+ // {\r
+ // double MinAb = min(AbA, AbB);\r
+ // double Ratio = MinAb/AbQ;\r
+ // double t = Ratio - opt_abx;\r
+ // // double Factor = 2.0/(1.0 + exp(-t));\r
+ // double Factor = min(Ratio, opt_abx)/opt_abx;\r
+ // if (opt_verbose)\r
+ // Log("Score %.4f Ab factor %.4f >%s\n", Hit.Score, Factor, QLabel.c_str());\r
+ // Hit.Score *= Factor;\r
+ // }\r
+ // }\r
+\r
+ extern FILE *g_fUChimeAlns;\r
+ if (g_fUChimeAlns != 0 && Hit.Div > 0.0)\r
+ {\r
+ void WriteChimeHitX(FILE *f, const ChimeHit2 &Hit);\r
+ WriteChimeHitX(g_fUChimeAlns, Hit);\r
+ }\r
+ }\r
+\r
+void AlignChime3(const string &Q3, const string &A3, const string &B3,\r
+ const string &QLabel, const string &ALabel, const string &BLabel,\r
+ ChimeHit2 &Hit)\r
+ {\r
+ if (opt_ucl)\r
+ AlignChimeLocal3(Q3, A3, B3, QLabel, ALabel, BLabel, Hit);\r
+ else\r
+ AlignChimeGlobal3(Q3, A3, B3, QLabel, ALabel, BLabel, Hit);\r
+ }\r
+\r
+static void StripGaps(const byte *Seq, unsigned L, string &s)\r
+ {\r
+ s.clear();\r
+ for (unsigned i = 0; i < L; ++i)\r
+ {\r
+ char c = Seq[i];\r
+ if (!isgap(c))\r
+ s.push_back(c);\r
+ }\r
+ }\r
+\r
+static void StripGapsAlloc(const SeqData &SDIn, SeqData &SDOut)\r
+ {\r
+ SDOut = SDIn;\r
+ byte *s = myalloc(byte, SDIn.L);\r
+ unsigned k = 0;\r
+ for (unsigned i = 0; i < SDIn.L; ++i)\r
+ {\r
+ char c = SDIn.Seq[i];\r
+ if (!isgap(c))\r
+ s[k++] = toupper(c);\r
+ }\r
+ SDOut.Seq = s;\r
+ SDOut.L = k;\r
+ }\r
+\r
+void AlignChime(const SeqData &QSD, const SeqData &ASD, const SeqData &BSD,\r
+ const string &PathQA, const string &PathQB, ChimeHit2 &Hit)\r
+ {\r
+ //if (opt_ucl)\r
+ // {\r
+ // AlignChimeLocal(QSD, ASD, BSD, PathQA, PathQB, Hit);\r
+ // return;\r
+ // }\r
+\r
+ string Q3;\r
+ string A3;\r
+ string B3;\r
+ Make3Way(QSD, ASD, BSD, PathQA, PathQB, Q3, A3, B3);\r
+\r
+ AlignChime3(Q3, A3, B3, QSD.Label, ASD.Label, BSD.Label, Hit);\r
+ }\r
+\r
+void AlignChime3SDRealign(const SeqData &QSD3, const SeqData &ASD3, const SeqData &BSD3,\r
+ ChimeHit2 &Hit)\r
+ {\r
+ SeqData QSD;\r
+ SeqData ASD;\r
+ SeqData BSD;\r
+ StripGapsAlloc(QSD3, QSD);\r
+ StripGapsAlloc(ASD3, ASD);\r
+ StripGapsAlloc(BSD3, BSD);\r
+\r
+ string PathQA;\r
+ string PathQB;\r
+ bool FoundQA = GlobalAlign(QSD, ASD, PathQA);\r
+ bool FoundQB = GlobalAlign(QSD, BSD, PathQB);\r
+ if (!FoundQA || !FoundQB)\r
+ {\r
+ Hit.Clear();\r
+ Hit.QLabel = QSD3.Label;\r
+ return;\r
+ }\r
+\r
+ AlignChime(QSD, ASD, BSD, PathQA, PathQB, Hit);\r
+\r
+ myfree((void *) QSD.Seq);\r
+ myfree((void *) ASD.Seq);\r
+ myfree((void *) BSD.Seq);\r
+ }\r
+\r
+void AlignChime3SD(const SeqData &QSD3, const SeqData &ASD3, const SeqData &BSD3,\r
+ ChimeHit2 &Hit)\r
+ {\r
+ if (opt_realign)\r
+ {\r
+ AlignChime3SDRealign(QSD3, ASD3, BSD3, Hit);\r
+ return;\r
+ }\r
+\r
+ string Q3;\r
+ string A3;\r
+ string B3;\r
+\r
+ const unsigned ColCount = QSD3.L;\r
+ asserta(ASD3.L == ColCount && BSD3.L == ColCount);\r
+\r
+ Q3.reserve(ColCount);\r
+ A3.reserve(ColCount);\r
+ B3.reserve(ColCount);\r
+\r
+ const byte *QS = QSD3.Seq;\r
+ const byte *AS = ASD3.Seq;\r
+ const byte *BS = BSD3.Seq;\r
+ for (unsigned Col = 0; Col < ColCount; ++Col)\r
+ {\r
+ byte q = toupper(QS[Col]);\r
+ byte a = toupper(AS[Col]);\r
+ byte b = toupper(BS[Col]);\r
+\r
+ if (isgap(q) && isgap(a) && isgap(b))\r
+ continue;\r
+\r
+ Q3.push_back(q);\r
+ A3.push_back(a);\r
+ B3.push_back(b);\r
+ }\r
+\r
+ AlignChime3(Q3, A3, B3, QSD3.Label, ASD3.Label, BSD3.Label, Hit);\r
+ }\r
projects / mothur.git / blobdiff