+++ /dev/null
-#include "muscle.h"\r
-#include "msa.h"\r
-#include "seqvect.h"\r
-#include "profile.h"\r
-#include "tree.h"\r
-\r
-// These global variables are a hack to allow the tree\r
-// dependent iteration code to communicate the edge\r
-// used to divide the tree. The three-way weighting\r
-// scheme needs to know this edge in order to compute\r
-// sequence weights.\r
-static const Tree *g_ptrMuscleTree = 0;\r
-unsigned g_uTreeSplitNode1 = NULL_NEIGHBOR;\r
-unsigned g_uTreeSplitNode2 = NULL_NEIGHBOR;\r
-\r
-void MSA::GetFractionalWeightedCounts(unsigned uColIndex, bool bNormalize,\r
- FCOUNT fcCounts[], FCOUNT *ptrfcGapStart, FCOUNT *ptrfcGapEnd,\r
- FCOUNT *ptrfcGapExtend, FCOUNT *ptrfOcc,\r
- FCOUNT *ptrfcLL, FCOUNT *ptrfcLG, FCOUNT *ptrfcGL, FCOUNT *ptrfcGG) const\r
- {\r
- const unsigned uSeqCount = GetSeqCount();\r
- const unsigned uColCount = GetColCount();\r
-\r
- memset(fcCounts, 0, g_AlphaSize*sizeof(FCOUNT));\r
- WEIGHT wTotal = 0;\r
- FCOUNT fGap = 0;\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- {\r
- const WEIGHT w = GetSeqWeight(uSeqIndex);\r
- if (IsGap(uSeqIndex, uColIndex))\r
- {\r
- fGap += w;\r
- continue;\r
- }\r
- else if (IsWildcard(uSeqIndex, uColIndex))\r
- {\r
- const unsigned uLetter = GetLetterEx(uSeqIndex, uColIndex);\r
- switch (g_Alpha)\r
- {\r
- case ALPHA_Amino:\r
- switch (uLetter)\r
- {\r
- case AX_B: // D or N\r
- fcCounts[AX_D] += w/2;\r
- fcCounts[AX_N] += w/2;\r
- break;\r
- case AX_Z: // E or Q\r
- fcCounts[AX_E] += w/2;\r
- fcCounts[AX_Q] += w/2;\r
- break;\r
- default: // any\r
- {\r
- const FCOUNT f = w/20;\r
- for (unsigned uLetter = 0; uLetter < 20; ++uLetter)\r
- fcCounts[uLetter] += f;\r
- break;\r
- }\r
- }\r
- break;\r
-\r
- case ALPHA_DNA:\r
- case ALPHA_RNA:\r
- switch (uLetter)\r
- {\r
- case AX_R: // G or A\r
- fcCounts[NX_G] += w/2;\r
- fcCounts[NX_A] += w/2;\r
- break;\r
- case AX_Y: // C or T/U\r
- fcCounts[NX_C] += w/2;\r
- fcCounts[NX_T] += w/2;\r
- break;\r
- default: // any\r
- const FCOUNT f = w/20;\r
- for (unsigned uLetter = 0; uLetter < 4; ++uLetter)\r
- fcCounts[uLetter] += f;\r
- break;\r
- }\r
- break;\r
-\r
- default:\r
- Quit("Alphabet %d not supported", g_Alpha);\r
- }\r
- continue;\r
- }\r
- unsigned uLetter = GetLetter(uSeqIndex, uColIndex);\r
- fcCounts[uLetter] += w;\r
- wTotal += w;\r
- }\r
- *ptrfOcc = (float) (1.0 - fGap);\r
-\r
- if (bNormalize && wTotal > 0)\r
- {\r
- if (wTotal > 1.001)\r
- Quit("wTotal=%g\n", wTotal);\r
- for (unsigned uLetter = 0; uLetter < g_AlphaSize; ++uLetter)\r
- fcCounts[uLetter] /= wTotal;\r
-// AssertNormalized(fcCounts);\r
- }\r
-\r
- FCOUNT fcStartCount = 0;\r
- if (uColIndex == 0)\r
- {\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- if (IsGap(uSeqIndex, uColIndex))\r
- fcStartCount += GetSeqWeight(uSeqIndex);\r
- }\r
- else\r
- {\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- if (IsGap(uSeqIndex, uColIndex) && !IsGap(uSeqIndex, uColIndex - 1))\r
- fcStartCount += GetSeqWeight(uSeqIndex);\r
- }\r
-\r
- FCOUNT fcEndCount = 0;\r
- if (uColCount - 1 == uColIndex)\r
- {\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- if (IsGap(uSeqIndex, uColIndex))\r
- fcEndCount += GetSeqWeight(uSeqIndex);\r
- }\r
- else\r
- {\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- if (IsGap(uSeqIndex, uColIndex) && !IsGap(uSeqIndex, uColIndex + 1))\r
- fcEndCount += GetSeqWeight(uSeqIndex);\r
- }\r
-\r
- FCOUNT LL = 0;\r
- FCOUNT LG = 0;\r
- FCOUNT GL = 0;\r
- FCOUNT GG = 0;\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- {\r
- WEIGHT w = GetSeqWeight(uSeqIndex);\r
- bool bLetterHere = !IsGap(uSeqIndex, uColIndex);\r
- bool bLetterPrev = (uColIndex == 0 || !IsGap(uSeqIndex, uColIndex - 1));\r
- if (bLetterHere)\r
- {\r
- if (bLetterPrev)\r
- LL += w;\r
- else\r
- GL += w;\r
- }\r
- else\r
- {\r
- if (bLetterPrev)\r
- LG += w;\r
- else\r
- GG += w;\r
- }\r
- }\r
-\r
- FCOUNT fcExtendCount = 0;\r
- if (uColIndex > 0 && uColIndex < GetColCount() - 1)\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- if (IsGap(uSeqIndex, uColIndex) && IsGap(uSeqIndex, uColIndex - 1) &&\r
- IsGap(uSeqIndex, uColIndex + 1))\r
- fcExtendCount += GetSeqWeight(uSeqIndex);\r
-\r
- *ptrfcLL = LL;\r
- *ptrfcLG = LG;\r
- *ptrfcGL = GL;\r
- *ptrfcGG = GG;\r
- *ptrfcGapStart = fcStartCount;\r
- *ptrfcGapEnd = fcEndCount;\r
- *ptrfcGapExtend = fcExtendCount;\r
- }\r
-\r
-// Return true if the given column has no gaps and all\r
-// its residues are in the same biochemical group.\r
-bool MSAColIsConservative(const MSA &msa, unsigned uColIndex)\r
- {\r
- extern unsigned ResidueGroup[];\r
-\r
- const unsigned uSeqCount = msa.GetColCount();\r
- if (0 == uSeqCount)\r
- Quit("MSAColIsConservative: empty alignment");\r
-\r
- if (msa.IsGap(0, uColIndex))\r
- return false;\r
-\r
- unsigned uLetter = msa.GetLetterEx(0, uColIndex);\r
- const unsigned uGroup = ResidueGroup[uLetter];\r
-\r
- for (unsigned uSeqIndex = 1; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- {\r
- if (msa.IsGap(uSeqIndex, uColIndex))\r
- return false;\r
- uLetter = msa.GetLetter(uSeqIndex, uColIndex);\r
- if (ResidueGroup[uLetter] != uGroup)\r
- return false;\r
- }\r
- return true;\r
- }\r
-\r
-void MSAFromSeqRange(const MSA &msaIn, unsigned uFromSeqIndex, unsigned uSeqCount,\r
- MSA &msaOut)\r
- {\r
- const unsigned uColCount = msaIn.GetColCount();\r
- msaOut.SetSize(uSeqCount, uColCount);\r
-\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- {\r
- const char *ptrName = msaIn.GetSeqName(uFromSeqIndex + uSeqIndex);\r
- msaOut.SetSeqName(uSeqIndex, ptrName);\r
-\r
- for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
- {\r
- const char c = msaIn.GetChar(uFromSeqIndex + uSeqIndex, uColIndex);\r
- msaOut.SetChar(uSeqIndex, uColIndex, c);\r
- }\r
- }\r
- }\r
-\r
-void MSAFromColRange(const MSA &msaIn, unsigned uFromColIndex, unsigned uColCount,\r
- MSA &msaOut)\r
- {\r
- const unsigned uSeqCount = msaIn.GetSeqCount();\r
- const unsigned uInColCount = msaIn.GetColCount();\r
-\r
- if (uFromColIndex + uColCount - 1 > uInColCount)\r
- Quit("MSAFromColRange, out of bounds");\r
-\r
- msaOut.SetSize(uSeqCount, uColCount);\r
-\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- {\r
- const char *ptrName = msaIn.GetSeqName(uSeqIndex);\r
- unsigned uId = msaIn.GetSeqId(uSeqIndex);\r
- msaOut.SetSeqName(uSeqIndex, ptrName);\r
- msaOut.SetSeqId(uSeqIndex, uId);\r
-\r
- for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
- {\r
- const char c = msaIn.GetChar(uSeqIndex, uFromColIndex + uColIndex);\r
- msaOut.SetChar(uSeqIndex, uColIndex, c);\r
- }\r
- }\r
- }\r
-\r
-void SeqVectFromMSA(const MSA &msa, SeqVect &v)\r
- {\r
- v.Clear();\r
- const unsigned uSeqCount = msa.GetSeqCount();\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- {\r
- Seq s;\r
- msa.GetSeq(uSeqIndex, s);\r
-\r
- s.StripGaps();\r
- //if (0 == s.Length())\r
- // continue;\r
-\r
- const char *ptrName = msa.GetSeqName(uSeqIndex);\r
- s.SetName(ptrName);\r
-\r
- v.AppendSeq(s);\r
- }\r
- }\r
-\r
-void DeleteGappedCols(MSA &msa)\r
- {\r
- unsigned uColIndex = 0;\r
- for (;;)\r
- {\r
- if (uColIndex >= msa.GetColCount())\r
- break;\r
- if (msa.IsGapColumn(uColIndex))\r
- msa.DeleteCol(uColIndex);\r
- else\r
- ++uColIndex;\r
- }\r
- }\r
-\r
-void MSAFromSeqSubset(const MSA &msaIn, const unsigned uSeqIndexes[], unsigned uSeqCount,\r
- MSA &msaOut)\r
- {\r
- const unsigned uColCount = msaIn.GetColCount();\r
- msaOut.SetSize(uSeqCount, uColCount);\r
- for (unsigned uSeqIndexOut = 0; uSeqIndexOut < uSeqCount; ++uSeqIndexOut)\r
- {\r
- unsigned uSeqIndexIn = uSeqIndexes[uSeqIndexOut];\r
- const char *ptrName = msaIn.GetSeqName(uSeqIndexIn);\r
- unsigned uId = msaIn.GetSeqId(uSeqIndexIn);\r
- msaOut.SetSeqName(uSeqIndexOut, ptrName);\r
- msaOut.SetSeqId(uSeqIndexOut, uId);\r
- for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
- {\r
- const char c = msaIn.GetChar(uSeqIndexIn, uColIndex);\r
- msaOut.SetChar(uSeqIndexOut, uColIndex, c);\r
- }\r
- }\r
- }\r
-\r
-void AssertMSAEqIgnoreCaseAndGaps(const MSA &msa1, const MSA &msa2)\r
- {\r
- const unsigned uSeqCount1 = msa1.GetSeqCount();\r
- const unsigned uSeqCount2 = msa2.GetSeqCount();\r
- if (uSeqCount1 != uSeqCount2)\r
- Quit("Seq count differs");\r
-\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount1; ++uSeqIndex)\r
- {\r
- Seq seq1;\r
- msa1.GetSeq(uSeqIndex, seq1);\r
-\r
- unsigned uId = msa1.GetSeqId(uSeqIndex);\r
- unsigned uSeqIndex2 = msa2.GetSeqIndex(uId);\r
-\r
- Seq seq2;\r
- msa2.GetSeq(uSeqIndex2, seq2);\r
-\r
- if (!seq1.EqIgnoreCaseAndGaps(seq2))\r
- {\r
- Log("Input:\n");\r
- seq1.LogMe();\r
- Log("Output:\n");\r
- seq2.LogMe();\r
- Quit("Seq %s differ ", msa1.GetSeqName(uSeqIndex));\r
- }\r
- }\r
- }\r
-\r
-void AssertMSAEq(const MSA &msa1, const MSA &msa2)\r
- {\r
- const unsigned uSeqCount1 = msa1.GetSeqCount();\r
- const unsigned uSeqCount2 = msa2.GetSeqCount();\r
- if (uSeqCount1 != uSeqCount2)\r
- Quit("Seq count differs");\r
-\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount1; ++uSeqIndex)\r
- {\r
- Seq seq1;\r
- msa1.GetSeq(uSeqIndex, seq1);\r
-\r
- unsigned uId = msa1.GetSeqId(uSeqIndex);\r
- unsigned uSeqIndex2 = msa2.GetSeqIndex(uId);\r
-\r
- Seq seq2;\r
- msa2.GetSeq(uSeqIndex2, seq2);\r
-\r
- if (!seq1.Eq(seq2))\r
- {\r
- Log("Input:\n");\r
- seq1.LogMe();\r
- Log("Output:\n");\r
- seq2.LogMe();\r
- Quit("Seq %s differ ", msa1.GetSeqName(uSeqIndex));\r
- }\r
- }\r
- }\r
-\r
-void SetMSAWeightsMuscle(MSA &msa)\r
- {\r
- SEQWEIGHT Method = GetSeqWeightMethod();\r
- switch (Method)\r
- {\r
- case SEQWEIGHT_None:\r
- msa.SetUniformWeights();\r
- return;\r
-\r
- case SEQWEIGHT_Henikoff:\r
- msa.SetHenikoffWeights();\r
- return;\r
-\r
- case SEQWEIGHT_HenikoffPB:\r
- msa.SetHenikoffWeightsPB();\r
- return;\r
-\r
- case SEQWEIGHT_GSC:\r
- msa.SetGSCWeights();\r
- return;\r
-\r
- case SEQWEIGHT_ClustalW:\r
- SetClustalWWeightsMuscle(msa);\r
- return;\r
- \r
- case SEQWEIGHT_ThreeWay:\r
- SetThreeWayWeightsMuscle(msa);\r
- return;\r
- }\r
- Quit("SetMSAWeightsMuscle, Invalid method=%d", Method);\r
- }\r
-\r
-static WEIGHT *g_MuscleWeights;\r
-static unsigned g_uMuscleIdCount;\r
-\r
-WEIGHT GetMuscleSeqWeightById(unsigned uId)\r
- {\r
- if (0 == g_MuscleWeights)\r
- Quit("g_MuscleWeights = 0");\r
- if (uId >= g_uMuscleIdCount)\r
- Quit("GetMuscleSeqWeightById(%u): count=%u",\r
- uId, g_uMuscleIdCount);\r
-\r
- return g_MuscleWeights[uId];\r
- }\r
-\r
-void SetMuscleTree(const Tree &tree)\r
- {\r
- g_ptrMuscleTree = &tree;\r
-\r
- if (SEQWEIGHT_ClustalW != GetSeqWeightMethod())\r
- return;\r
-\r
- delete[] g_MuscleWeights;\r
-\r
- const unsigned uLeafCount = tree.GetLeafCount();\r
- g_uMuscleIdCount = uLeafCount;\r
- g_MuscleWeights = new WEIGHT[uLeafCount];\r
- CalcClustalWWeights(tree, g_MuscleWeights);\r
- }\r
-\r
-void SetClustalWWeightsMuscle(MSA &msa)\r
- {\r
- if (0 == g_MuscleWeights)\r
- Quit("g_MuscleWeights = 0");\r
- const unsigned uSeqCount = msa.GetSeqCount();\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- {\r
- const unsigned uId = msa.GetSeqId(uSeqIndex);\r
- if (uId >= g_uMuscleIdCount)\r
- Quit("SetClustalWWeightsMuscle: id out of range");\r
- msa.SetSeqWeight(uSeqIndex, g_MuscleWeights[uId]);\r
- }\r
- msa.NormalizeWeights((WEIGHT) 1.0);\r
- }\r
-\r
-#define LOCAL_VERBOSE 0\r
-\r
-void SetThreeWayWeightsMuscle(MSA &msa)\r
- {\r
- if (NULL_NEIGHBOR == g_uTreeSplitNode1 || NULL_NEIGHBOR == g_uTreeSplitNode2)\r
- {\r
- msa.SetHenikoffWeightsPB();\r
- return;\r
- }\r
-\r
- const unsigned uMuscleSeqCount = g_ptrMuscleTree->GetLeafCount();\r
- WEIGHT *Weights = new WEIGHT[uMuscleSeqCount];\r
-\r
- CalcThreeWayWeights(*g_ptrMuscleTree, g_uTreeSplitNode1, g_uTreeSplitNode2,\r
- Weights);\r
-\r
- const unsigned uMSASeqCount = msa.GetSeqCount();\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uMSASeqCount; ++uSeqIndex)\r
- {\r
- const unsigned uId = msa.GetSeqId(uSeqIndex);\r
- if (uId >= uMuscleSeqCount)\r
- Quit("SetThreeWayWeightsMuscle: id out of range");\r
- msa.SetSeqWeight(uSeqIndex, Weights[uId]);\r
- }\r
-#if LOCAL_VERBOSE\r
- {\r
- Log("SetThreeWayWeightsMuscle\n");\r
- for (unsigned n = 0; n < uMSASeqCount; ++n)\r
- {\r
- const unsigned uId = msa.GetSeqId(n);\r
- Log("%20.20s %6.3f\n", msa.GetSeqName(n), Weights[uId]);\r
- }\r
- }\r
-#endif\r
- msa.NormalizeWeights((WEIGHT) 1.0);\r
-\r
- delete[] Weights;\r
- }\r
-\r
-// Append msa2 at the end of msa1\r
-void MSAAppend(MSA &msa1, const MSA &msa2)\r
- {\r
- const unsigned uSeqCount = msa1.GetSeqCount();\r
-\r
- const unsigned uColCount1 = msa1.GetColCount();\r
- const unsigned uColCount2 = msa2.GetColCount();\r
- const unsigned uColCountCat = uColCount1 + uColCount2;\r
-\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- {\r
- unsigned uId = msa1.GetSeqId(uSeqIndex);\r
- unsigned uSeqIndex2 = msa2.GetSeqIndex(uId);\r
- for (unsigned uColIndex = 0; uColIndex < uColCount2; ++uColIndex)\r
- {\r
- const char c = msa2.GetChar(uSeqIndex2, uColIndex);\r
- msa1.SetChar(uSeqIndex, uColCount1 + uColIndex, c);\r
- }\r
- }\r
- }\r
-\r
-// "Catenate" two MSAs (by bad analogy with UNIX cat command).\r
-// msa1 and msa2 must have same sequence names, but possibly\r
-// in a different order.\r
-// msaCat is the combined alignment produce by appending\r
-// sequences in msa2 to sequences in msa1.\r
-void MSACat(const MSA &msa1, const MSA &msa2, MSA &msaCat)\r
- {\r
- const unsigned uSeqCount = msa1.GetSeqCount();\r
-\r
- const unsigned uColCount1 = msa1.GetColCount();\r
- const unsigned uColCount2 = msa2.GetColCount();\r
- const unsigned uColCountCat = uColCount1 + uColCount2;\r
-\r
- msaCat.SetSize(uSeqCount, uColCountCat);\r
-\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- {\r
- for (unsigned uColIndex = 0; uColIndex < uColCount1; ++uColIndex)\r
- {\r
- const char c = msa1.GetChar(uSeqIndex, uColIndex);\r
- msaCat.SetChar(uSeqIndex, uColIndex, c);\r
- }\r
-\r
- const char *ptrSeqName = msa1.GetSeqName(uSeqIndex);\r
- unsigned uSeqIndex2;\r
- msaCat.SetSeqName(uSeqIndex, ptrSeqName);\r
- bool bFound = msa2.GetSeqIndex(ptrSeqName, &uSeqIndex2);\r
- if (bFound)\r
- {\r
- for (unsigned uColIndex = 0; uColIndex < uColCount2; ++uColIndex)\r
- {\r
- const char c = msa2.GetChar(uSeqIndex2, uColIndex);\r
- msaCat.SetChar(uSeqIndex, uColCount1 + uColIndex, c);\r
- }\r
- }\r
- else\r
- {\r
- for (unsigned uColIndex = 0; uColIndex < uColCount2; ++uColIndex)\r
- msaCat.SetChar(uSeqIndex, uColCount1 + uColIndex, '-');\r
- }\r
- }\r
- }\r
git://source.jalview.org / jabaws.git / blobdiff