+++ /dev/null
-#include "muscle.h"\r
-#include "msa.h"\r
-#include "textfile.h"\r
-#include "seq.h"\r
-#include <math.h>\r
-\r
-const unsigned DEFAULT_SEQ_LENGTH = 500;\r
-\r
-unsigned MSA::m_uIdCount = 0;\r
-\r
-MSA::MSA()\r
- {\r
- m_uSeqCount = 0;\r
- m_uColCount = 0;\r
-\r
- m_szSeqs = 0;\r
- m_szNames = 0;\r
- m_Weights = 0;\r
-\r
- m_IdToSeqIndex = 0;\r
- m_SeqIndexToId = 0;\r
-\r
- m_uCacheSeqCount = 0;\r
- m_uCacheSeqLength = 0;\r
- }\r
-\r
-MSA::~MSA()\r
- {\r
- Free();\r
- }\r
-\r
-void MSA::Free()\r
- {\r
- for (unsigned n = 0; n < m_uSeqCount; ++n)\r
- {\r
- delete[] m_szSeqs[n];\r
- delete[] m_szNames[n];\r
- }\r
-\r
- delete[] m_szSeqs;\r
- delete[] m_szNames;\r
- delete[] m_Weights;\r
- delete[] m_IdToSeqIndex;\r
- delete[] m_SeqIndexToId;\r
-\r
- m_uSeqCount = 0;\r
- m_uColCount = 0;\r
-\r
- m_szSeqs = 0;\r
- m_szNames = 0;\r
- m_Weights = 0;\r
-\r
- m_IdToSeqIndex = 0;\r
- m_SeqIndexToId = 0;\r
- }\r
-\r
-void MSA::SetSize(unsigned uSeqCount, unsigned uColCount)\r
- {\r
- Free();\r
-\r
- m_uSeqCount = uSeqCount;\r
- m_uCacheSeqLength = uColCount;\r
- m_uColCount = 0;\r
-\r
- if (0 == uSeqCount && 0 == uColCount)\r
- return;\r
-\r
- m_szSeqs = new char *[uSeqCount];\r
- m_szNames = new char *[uSeqCount];\r
- m_Weights = new WEIGHT[uSeqCount];\r
-\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- {\r
- m_szSeqs[uSeqIndex] = new char[uColCount+1];\r
- m_szNames[uSeqIndex] = 0;\r
-#if DEBUG\r
- m_Weights[uSeqIndex] = BTInsane;\r
- memset(m_szSeqs[uSeqIndex], '?', uColCount);\r
-#endif\r
- m_szSeqs[uSeqIndex][uColCount] = 0;\r
- }\r
-\r
- if (m_uIdCount > 0)\r
- {\r
- m_IdToSeqIndex = new unsigned[m_uIdCount];\r
- m_SeqIndexToId = new unsigned[m_uSeqCount];\r
-#if DEBUG\r
- memset(m_IdToSeqIndex, 0xff, m_uIdCount*sizeof(unsigned));\r
- memset(m_SeqIndexToId, 0xff, m_uSeqCount*sizeof(unsigned));\r
-#endif\r
- }\r
- }\r
-\r
-void MSA::LogMe() const\r
- {\r
- if (0 == GetColCount())\r
- {\r
- Log("MSA empty\n");\r
- return;\r
- }\r
-\r
- const unsigned uColsPerLine = 50;\r
- unsigned uLinesPerSeq = (GetColCount() - 1)/uColsPerLine + 1;\r
- for (unsigned n = 0; n < uLinesPerSeq; ++n)\r
- {\r
- unsigned i;\r
- unsigned iStart = n*uColsPerLine;\r
- unsigned iEnd = GetColCount();\r
- if (iEnd - iStart + 1 > uColsPerLine)\r
- iEnd = iStart + uColsPerLine;\r
- Log(" ");\r
- for (i = iStart; i < iEnd; ++i)\r
- Log("%u", i%10);\r
- Log("\n");\r
- Log(" ");\r
- for (i = iStart; i + 9 < iEnd; i += 10)\r
- Log("%-10u", i);\r
- if (n == uLinesPerSeq - 1)\r
- Log(" %-10u", GetColCount());\r
- Log("\n");\r
- for (unsigned uSeqIndex = 0; uSeqIndex < m_uSeqCount; ++uSeqIndex)\r
- {\r
- Log("%12.12s", m_szNames[uSeqIndex]);\r
- if (m_Weights[uSeqIndex] != BTInsane)\r
- Log(" (%5.3f)", m_Weights[uSeqIndex]);\r
- else\r
- Log(" ");\r
- Log(" ");\r
- for (i = iStart; i < iEnd; ++i)\r
- Log("%c", GetChar(uSeqIndex, i));\r
- if (0 != m_SeqIndexToId)\r
- Log(" [%5u]", m_SeqIndexToId[uSeqIndex]);\r
- Log("\n");\r
- }\r
- Log("\n\n");\r
- }\r
- }\r
-\r
-char MSA::GetChar(unsigned uSeqIndex, unsigned uIndex) const\r
- {\r
-// TODO: Performance cost?\r
- if (uSeqIndex >= m_uSeqCount || uIndex >= m_uColCount)\r
- Quit("MSA::GetChar(%u/%u,%u/%u)",\r
- uSeqIndex, m_uSeqCount, uIndex, m_uColCount);\r
-\r
- char c = m_szSeqs[uSeqIndex][uIndex];\r
-// assert(IsLegalChar(c));\r
- return c;\r
- }\r
-\r
-unsigned MSA::GetLetter(unsigned uSeqIndex, unsigned uIndex) const\r
- {\r
-// TODO: Performance cost?\r
- char c = GetChar(uSeqIndex, uIndex);\r
- unsigned uLetter = CharToLetter(c);\r
- if (uLetter >= 20)\r
- {\r
- char c = ' ';\r
- if (uSeqIndex < m_uSeqCount && uIndex < m_uColCount)\r
- c = m_szSeqs[uSeqIndex][uIndex];\r
- Quit("MSA::GetLetter(%u/%u, %u/%u)='%c'/%u",\r
- uSeqIndex, m_uSeqCount, uIndex, m_uColCount, c, uLetter);\r
- }\r
- return uLetter;\r
- }\r
-\r
-unsigned MSA::GetLetterEx(unsigned uSeqIndex, unsigned uIndex) const\r
- {\r
-// TODO: Performance cost?\r
- char c = GetChar(uSeqIndex, uIndex);\r
- unsigned uLetter = CharToLetterEx(c);\r
- return uLetter;\r
- }\r
-\r
-void MSA::SetSeqName(unsigned uSeqIndex, const char szName[])\r
- {\r
- if (uSeqIndex >= m_uSeqCount)\r
- Quit("MSA::SetSeqName(%u, %s), count=%u", uSeqIndex, m_uSeqCount);\r
- delete[] m_szNames[uSeqIndex];\r
- int n = (int) strlen(szName) + 1;\r
- m_szNames[uSeqIndex] = new char[n];\r
- memcpy(m_szNames[uSeqIndex], szName, n);\r
- }\r
-\r
-const char *MSA::GetSeqName(unsigned uSeqIndex) const\r
- {\r
- if (uSeqIndex >= m_uSeqCount)\r
- Quit("MSA::GetSeqName(%u), count=%u", uSeqIndex, m_uSeqCount);\r
- return m_szNames[uSeqIndex];\r
- }\r
-\r
-bool MSA::IsGap(unsigned uSeqIndex, unsigned uIndex) const\r
- {\r
- char c = GetChar(uSeqIndex, uIndex);\r
- return IsGapChar(c);\r
- }\r
-\r
-bool MSA::IsWildcard(unsigned uSeqIndex, unsigned uIndex) const\r
- {\r
- char c = GetChar(uSeqIndex, uIndex);\r
- return IsWildcardChar(c);\r
- }\r
-\r
-void MSA::SetChar(unsigned uSeqIndex, unsigned uIndex, char c)\r
- {\r
- if (uSeqIndex >= m_uSeqCount || uIndex > m_uCacheSeqLength)\r
- Quit("MSA::SetChar(%u,%u)", uSeqIndex, uIndex);\r
-\r
- if (uIndex == m_uCacheSeqLength)\r
- {\r
- const unsigned uNewCacheSeqLength = m_uCacheSeqLength + DEFAULT_SEQ_LENGTH;\r
- for (unsigned n = 0; n < m_uSeqCount; ++n)\r
- {\r
- char *ptrNewSeq = new char[uNewCacheSeqLength+1];\r
- memcpy(ptrNewSeq, m_szSeqs[n], m_uCacheSeqLength);\r
- memset(ptrNewSeq + m_uCacheSeqLength, '?', DEFAULT_SEQ_LENGTH);\r
- ptrNewSeq[uNewCacheSeqLength] = 0;\r
- delete[] m_szSeqs[n];\r
- m_szSeqs[n] = ptrNewSeq;\r
- }\r
-\r
- m_uColCount = uIndex;\r
- m_uCacheSeqLength = uNewCacheSeqLength;\r
- }\r
-\r
- if (uIndex >= m_uColCount)\r
- m_uColCount = uIndex + 1;\r
- m_szSeqs[uSeqIndex][uIndex] = c;\r
- }\r
-\r
-void MSA::GetSeq(unsigned uSeqIndex, Seq &seq) const\r
- {\r
- assert(uSeqIndex < m_uSeqCount);\r
-\r
- seq.Clear();\r
-\r
- for (unsigned n = 0; n < m_uColCount; ++n)\r
- if (!IsGap(uSeqIndex, n))\r
- {\r
- char c = GetChar(uSeqIndex, n);\r
- if (!isalpha(c))\r
- Quit("Invalid character '%c' in sequence", c);\r
- c = toupper(c);\r
- seq.push_back(c);\r
- }\r
- const char *ptrName = GetSeqName(uSeqIndex);\r
- seq.SetName(ptrName);\r
- }\r
-\r
-bool MSA::HasGap() const\r
- {\r
- for (unsigned uSeqIndex = 0; uSeqIndex < GetSeqCount(); ++uSeqIndex)\r
- for (unsigned n = 0; n < GetColCount(); ++n)\r
- if (IsGap(uSeqIndex, n))\r
- return true;\r
- return false;\r
- }\r
-\r
-bool MSA::IsLegalLetter(unsigned uLetter) const\r
- {\r
- return uLetter < 20;\r
- }\r
-\r
-void MSA::SetSeqCount(unsigned uSeqCount)\r
- {\r
- Free();\r
- SetSize(uSeqCount, DEFAULT_SEQ_LENGTH);\r
- }\r
-\r
-void MSA::CopyCol(unsigned uFromCol, unsigned uToCol)\r
- {\r
- assert(uFromCol < GetColCount());\r
- assert(uToCol < GetColCount());\r
- if (uFromCol == uToCol)\r
- return;\r
-\r
- for (unsigned uSeqIndex = 0; uSeqIndex < GetSeqCount(); ++uSeqIndex)\r
- {\r
- const char c = GetChar(uSeqIndex, uFromCol);\r
- SetChar(uSeqIndex, uToCol, c);\r
- }\r
- }\r
-\r
-void MSA::Copy(const MSA &msa)\r
- {\r
- Free();\r
- const unsigned uSeqCount = msa.GetSeqCount();\r
- const unsigned uColCount = msa.GetColCount();\r
- SetSize(uSeqCount, uColCount);\r
-\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- {\r
- SetSeqName(uSeqIndex, msa.GetSeqName(uSeqIndex));\r
- const unsigned uId = msa.GetSeqId(uSeqIndex);\r
- SetSeqId(uSeqIndex, uId);\r
- for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
- {\r
- const char c = msa.GetChar(uSeqIndex, uColIndex);\r
- SetChar(uSeqIndex, uColIndex, c);\r
- }\r
- }\r
- }\r
-\r
-bool MSA::IsGapColumn(unsigned uColIndex) const\r
- {\r
- assert(GetSeqCount() > 0);\r
- for (unsigned uSeqIndex = 0; uSeqIndex < GetSeqCount(); ++uSeqIndex)\r
- if (!IsGap(uSeqIndex, uColIndex))\r
- return false;\r
- return true;\r
- }\r
-\r
-bool MSA::GetSeqIndex(const char *ptrSeqName, unsigned *ptruSeqIndex) const\r
- {\r
- for (unsigned uSeqIndex = 0; uSeqIndex < GetSeqCount(); ++uSeqIndex)\r
- if (0 == stricmp(ptrSeqName, GetSeqName(uSeqIndex)))\r
- {\r
- *ptruSeqIndex = uSeqIndex;\r
- return true;\r
- }\r
- return false;\r
- }\r
-\r
-void MSA::DeleteCol(unsigned uColIndex)\r
- {\r
- assert(uColIndex < m_uColCount);\r
- size_t n = m_uColCount - uColIndex;\r
- if (n > 0)\r
- {\r
- for (unsigned uSeqIndex = 0; uSeqIndex < GetSeqCount(); ++uSeqIndex)\r
- {\r
- char *ptrSeq = m_szSeqs[uSeqIndex];\r
- memmove(ptrSeq + uColIndex, ptrSeq + uColIndex + 1, n);\r
- }\r
- }\r
- --m_uColCount;\r
- }\r
-\r
-void MSA::DeleteColumns(unsigned uColIndex, unsigned uColCount)\r
- {\r
- for (unsigned n = 0; n < uColCount; ++n)\r
- DeleteCol(uColIndex);\r
- }\r
-\r
-void MSA::FromFile(TextFile &File)\r
- {\r
- FromFASTAFile(File);\r
- }\r
-\r
-// Weights sum to 1, WCounts sum to NIC\r
-WEIGHT MSA::GetSeqWeight(unsigned uSeqIndex) const\r
- {\r
- assert(uSeqIndex < m_uSeqCount);\r
- WEIGHT w = m_Weights[uSeqIndex];\r
- if (w == wInsane)\r
- Quit("Seq weight not set");\r
- return w;\r
- }\r
-\r
-void MSA::SetSeqWeight(unsigned uSeqIndex, WEIGHT w) const\r
- {\r
- assert(uSeqIndex < m_uSeqCount);\r
- m_Weights[uSeqIndex] = w;\r
- }\r
-\r
-void MSA::NormalizeWeights(WEIGHT wDesiredTotal) const\r
- {\r
- WEIGHT wTotal = 0;\r
- for (unsigned uSeqIndex = 0; uSeqIndex < m_uSeqCount; ++uSeqIndex)\r
- wTotal += m_Weights[uSeqIndex];\r
-\r
- if (0 == wTotal)\r
- return;\r
-\r
- const WEIGHT f = wDesiredTotal/wTotal;\r
- for (unsigned uSeqIndex = 0; uSeqIndex < m_uSeqCount; ++uSeqIndex)\r
- m_Weights[uSeqIndex] *= f;\r
- }\r
-\r
-void MSA::CalcWeights() const\r
- {\r
- Quit("Calc weights not implemented");\r
- }\r
-\r
-static void FmtChar(char c, unsigned uWidth)\r
- {\r
- Log("%c", c);\r
- for (unsigned n = 0; n < uWidth - 1; ++n)\r
- Log(" ");\r
- }\r
-\r
-static void FmtInt(unsigned u, unsigned uWidth)\r
- {\r
- static char szStr[1024];\r
- assert(uWidth < sizeof(szStr));\r
- if (u > 0)\r
- sprintf(szStr, "%u", u);\r
- else\r
- strcpy(szStr, ".");\r
- Log(szStr);\r
- unsigned n = (unsigned) strlen(szStr);\r
- if (n < uWidth)\r
- for (unsigned i = 0; i < uWidth - n; ++i)\r
- Log(" ");\r
- }\r
-\r
-static void FmtInt0(unsigned u, unsigned uWidth)\r
- {\r
- static char szStr[1024];\r
- assert(uWidth < sizeof(szStr));\r
- sprintf(szStr, "%u", u);\r
- Log(szStr);\r
- unsigned n = (unsigned) strlen(szStr);\r
- if (n < uWidth)\r
- for (unsigned i = 0; i < uWidth - n; ++i)\r
- Log(" ");\r
- }\r
-\r
-static void FmtPad(unsigned n)\r
- {\r
- for (unsigned i = 0; i < n; ++i)\r
- Log(" ");\r
- }\r
-\r
-void MSA::FromSeq(const Seq &s)\r
- {\r
- unsigned uSeqLength = s.Length();\r
- SetSize(1, uSeqLength);\r
- SetSeqName(0, s.GetName());\r
- if (0 != m_SeqIndexToId)\r
- SetSeqId(0, s.GetId());\r
- for (unsigned n = 0; n < uSeqLength; ++n)\r
- SetChar(0, n, s[n]);\r
- }\r
-\r
-unsigned MSA::GetCharCount(unsigned uSeqIndex, unsigned uColIndex) const\r
- {\r
- assert(uSeqIndex < GetSeqCount());\r
- assert(uColIndex < GetColCount());\r
-\r
- unsigned uCol = 0;\r
- for (unsigned n = 0; n <= uColIndex; ++n)\r
- if (!IsGap(uSeqIndex, n))\r
- ++uCol;\r
- return uCol;\r
- }\r
-\r
-void MSA::CopySeq(unsigned uToSeqIndex, const MSA &msaFrom, unsigned uFromSeqIndex)\r
- {\r
- assert(uToSeqIndex < m_uSeqCount);\r
- const unsigned uColCount = msaFrom.GetColCount();\r
- assert(m_uColCount == uColCount ||\r
- (0 == m_uColCount && uColCount <= m_uCacheSeqLength));\r
-\r
- memcpy(m_szSeqs[uToSeqIndex], msaFrom.GetSeqBuffer(uFromSeqIndex), uColCount);\r
- SetSeqName(uToSeqIndex, msaFrom.GetSeqName(uFromSeqIndex));\r
- if (0 == m_uColCount)\r
- m_uColCount = uColCount;\r
- }\r
-\r
-const char *MSA::GetSeqBuffer(unsigned uSeqIndex) const\r
- {\r
- assert(uSeqIndex < m_uSeqCount);\r
- return m_szSeqs[uSeqIndex];\r
- }\r
-\r
-void MSA::DeleteSeq(unsigned uSeqIndex)\r
- {\r
- assert(uSeqIndex < m_uSeqCount);\r
-\r
- delete m_szSeqs[uSeqIndex];\r
- delete m_szNames[uSeqIndex];\r
-\r
- const unsigned uBytesToMove = (m_uSeqCount - uSeqIndex)*sizeof(char *);\r
- if (uBytesToMove > 0)\r
- {\r
- memmove(m_szSeqs + uSeqIndex, m_szSeqs + uSeqIndex + 1, uBytesToMove);\r
- memmove(m_szNames + uSeqIndex, m_szNames + uSeqIndex + 1, uBytesToMove);\r
- }\r
-\r
- --m_uSeqCount;\r
-\r
- delete[] m_Weights;\r
- m_Weights = 0;\r
- }\r
-\r
-bool MSA::IsEmptyCol(unsigned uColIndex) const\r
- {\r
- const unsigned uSeqCount = GetSeqCount();\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- if (!IsGap(uSeqIndex, uColIndex))\r
- return false;\r
- return true;\r
- }\r
-\r
-//void MSA::DeleteEmptyCols(bool bProgress)\r
-// {\r
-// unsigned uColCount = GetColCount();\r
-// for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
-// {\r
-// if (IsEmptyCol(uColIndex))\r
-// {\r
-// if (bProgress)\r
-// {\r
-// Log("Deleting col %u of %u\n", uColIndex, uColCount);\r
-// printf("Deleting col %u of %u\n", uColIndex, uColCount);\r
-// }\r
-// DeleteCol(uColIndex);\r
-// --uColCount;\r
-// }\r
-// }\r
-// }\r
-\r
-unsigned MSA::AlignedColIndexToColIndex(unsigned uAlignedColIndex) const\r
- {\r
- Quit("MSA::AlignedColIndexToColIndex not implemented");\r
- return 0;\r
- }\r
-\r
-WEIGHT MSA::GetTotalSeqWeight() const\r
- {\r
- WEIGHT wTotal = 0;\r
- const unsigned uSeqCount = GetSeqCount();\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- wTotal += m_Weights[uSeqIndex];\r
- return wTotal;\r
- }\r
-\r
-bool MSA::SeqsEq(const MSA &a1, unsigned uSeqIndex1, const MSA &a2,\r
- unsigned uSeqIndex2)\r
- {\r
- Seq s1;\r
- Seq s2;\r
-\r
- a1.GetSeq(uSeqIndex1, s1);\r
- a2.GetSeq(uSeqIndex2, s2);\r
-\r
- s1.StripGaps();\r
- s2.StripGaps();\r
-\r
- return s1.EqIgnoreCase(s2);\r
- }\r
-\r
-unsigned MSA::GetSeqLength(unsigned uSeqIndex) const\r
- {\r
- assert(uSeqIndex < GetSeqCount());\r
-\r
- const unsigned uColCount = GetColCount();\r
- unsigned uLength = 0;\r
- for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
- if (!IsGap(uSeqIndex, uColIndex))\r
- ++uLength;\r
- return uLength;\r
- }\r
-\r
-void MSA::GetPWID(unsigned uSeqIndex1, unsigned uSeqIndex2, double *ptrPWID,\r
- unsigned *ptruPosCount) const\r
- {\r
- assert(uSeqIndex1 < GetSeqCount());\r
- assert(uSeqIndex2 < GetSeqCount());\r
-\r
- unsigned uSameCount = 0;\r
- unsigned uPosCount = 0;\r
- const unsigned uColCount = GetColCount();\r
- for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
- {\r
- char c1 = GetChar(uSeqIndex1, uColIndex);\r
- if (IsGapChar(c1))\r
- continue;\r
- char c2 = GetChar(uSeqIndex2, uColIndex);\r
- if (IsGapChar(c2))\r
- continue;\r
- ++uPosCount;\r
- if (c1 == c2)\r
- ++uSameCount;\r
- }\r
- *ptruPosCount = uPosCount;\r
- if (uPosCount > 0)\r
- *ptrPWID = 100.0 * (double) uSameCount / (double) uPosCount;\r
- else\r
- *ptrPWID = 0;\r
- }\r
-\r
-void MSA::UnWeight()\r
- {\r
- for (unsigned uSeqIndex = 0; uSeqIndex < GetSeqCount(); ++uSeqIndex)\r
- m_Weights[uSeqIndex] = BTInsane;\r
- }\r
-\r
-unsigned MSA::UniqueResidueTypes(unsigned uColIndex) const\r
- {\r
- assert(uColIndex < GetColCount());\r
-\r
- unsigned Counts[MAX_ALPHA];\r
- memset(Counts, 0, sizeof(Counts));\r
- const unsigned uSeqCount = GetSeqCount();\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- {\r
- if (IsGap(uSeqIndex, uColIndex) || IsWildcard(uSeqIndex, uColIndex))\r
- continue;\r
- const unsigned uLetter = GetLetter(uSeqIndex, uColIndex);\r
- ++(Counts[uLetter]);\r
- }\r
- unsigned uUniqueCount = 0;\r
- for (unsigned uLetter = 0; uLetter < g_AlphaSize; ++uLetter)\r
- if (Counts[uLetter] > 0)\r
- ++uUniqueCount;\r
- return uUniqueCount;\r
- }\r
-\r
-double MSA::GetOcc(unsigned uColIndex) const\r
- {\r
- unsigned uGapCount = 0;\r
- for (unsigned uSeqIndex = 0; uSeqIndex < GetSeqCount(); ++uSeqIndex)\r
- if (IsGap(uSeqIndex, uColIndex))\r
- ++uGapCount;\r
- unsigned uSeqCount = GetSeqCount();\r
- return (double) (uSeqCount - uGapCount) / (double) uSeqCount;\r
- }\r
-\r
-void MSA::ToFile(TextFile &File) const\r
- {\r
- if (g_bMSF)\r
- ToMSFFile(File);\r
- else if (g_bAln)\r
- ToAlnFile(File);\r
- else if (g_bHTML)\r
- ToHTMLFile(File);\r
- else if (g_bPHYS)\r
- ToPhySequentialFile(File);\r
- else if (g_bPHYI)\r
- ToPhyInterleavedFile(File);\r
- else\r
- ToFASTAFile(File);\r
- if (0 != g_pstrScoreFileName)\r
- WriteScoreFile(*this);\r
- }\r
-\r
-bool MSA::ColumnHasGap(unsigned uColIndex) const\r
- {\r
- const unsigned uSeqCount = GetSeqCount();\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- if (IsGap(uSeqIndex, uColIndex))\r
- return true;\r
- return false;\r
- }\r
-\r
-void MSA::SetIdCount(unsigned uIdCount)\r
- {\r
- //if (m_uIdCount != 0)\r
- // Quit("MSA::SetIdCount: may only be called once");\r
-\r
- if (m_uIdCount > 0)\r
- {\r
- if (uIdCount > m_uIdCount)\r
- Quit("MSA::SetIdCount: cannot increase count");\r
- return;\r
- }\r
- m_uIdCount = uIdCount;\r
- }\r
-\r
-void MSA::SetSeqId(unsigned uSeqIndex, unsigned uId)\r
- {\r
- assert(uSeqIndex < m_uSeqCount);\r
- assert(uId < m_uIdCount);\r
- if (0 == m_SeqIndexToId)\r
- {\r
- if (0 == m_uIdCount)\r
- Quit("MSA::SetSeqId, SetIdCount has not been called");\r
- m_IdToSeqIndex = new unsigned[m_uIdCount];\r
- m_SeqIndexToId = new unsigned[m_uSeqCount];\r
-\r
- memset(m_IdToSeqIndex, 0xff, m_uIdCount*sizeof(unsigned));\r
- memset(m_SeqIndexToId, 0xff, m_uSeqCount*sizeof(unsigned));\r
- }\r
- m_SeqIndexToId[uSeqIndex] = uId;\r
- m_IdToSeqIndex[uId] = uSeqIndex;\r
- }\r
-\r
-unsigned MSA::GetSeqIndex(unsigned uId) const\r
- {\r
- assert(uId < m_uIdCount);\r
- assert(0 != m_IdToSeqIndex);\r
- unsigned uSeqIndex = m_IdToSeqIndex[uId];\r
- assert(uSeqIndex < m_uSeqCount);\r
- return uSeqIndex;\r
- }\r
-\r
-bool MSA::GetSeqIndex(unsigned uId, unsigned *ptruIndex) const\r
- {\r
- for (unsigned uSeqIndex = 0; uSeqIndex < m_uSeqCount; ++uSeqIndex)\r
- {\r
- if (uId == m_SeqIndexToId[uSeqIndex])\r
- {\r
- *ptruIndex = uSeqIndex;\r
- return true;\r
- }\r
- }\r
- return false;\r
- }\r
-\r
-unsigned MSA::GetSeqId(unsigned uSeqIndex) const\r
- {\r
- assert(uSeqIndex < m_uSeqCount);\r
- unsigned uId = m_SeqIndexToId[uSeqIndex];\r
- assert(uId < m_uIdCount);\r
- return uId;\r
- }\r
-\r
-bool MSA::WeightsSet() const\r
- {\r
- return BTInsane != m_Weights[0];\r
- }\r
-\r
-void MSASubsetByIds(const MSA &msaIn, const unsigned Ids[], unsigned uIdCount,\r
- MSA &msaOut)\r
- {\r
- const unsigned uColCount = msaIn.GetColCount();\r
- msaOut.SetSize(uIdCount, uColCount);\r
- for (unsigned uSeqIndexOut = 0; uSeqIndexOut < uIdCount; ++uSeqIndexOut)\r
- {\r
- const unsigned uId = Ids[uSeqIndexOut];\r
-\r
- const unsigned uSeqIndexIn = msaIn.GetSeqIndex(uId);\r
- const char *ptrName = msaIn.GetSeqName(uSeqIndexIn);\r
-\r
- msaOut.SetSeqId(uSeqIndexOut, uId);\r
- msaOut.SetSeqName(uSeqIndexOut, ptrName);\r
-\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
-// Caller must allocate ptrSeq and ptrLabel as new char[n].\r
-void MSA::AppendSeq(char *ptrSeq, unsigned uSeqLength, char *ptrLabel)\r
- {\r
- if (m_uSeqCount > m_uCacheSeqCount)\r
- Quit("Internal error MSA::AppendSeq");\r
- if (m_uSeqCount == m_uCacheSeqCount)\r
- ExpandCache(m_uSeqCount + 4, uSeqLength);\r
- m_szSeqs[m_uSeqCount] = ptrSeq;\r
- m_szNames[m_uSeqCount] = ptrLabel;\r
- ++m_uSeqCount;\r
- }\r
-\r
-void MSA::ExpandCache(unsigned uSeqCount, unsigned uColCount)\r
- {\r
- if (m_IdToSeqIndex != 0 || m_SeqIndexToId != 0 || uSeqCount < m_uSeqCount)\r
- Quit("Internal error MSA::ExpandCache");\r
-\r
- if (m_uSeqCount > 0 && uColCount != m_uColCount)\r
- Quit("Internal error MSA::ExpandCache, ColCount changed");\r
-\r
- char **NewSeqs = new char *[uSeqCount];\r
- char **NewNames = new char *[uSeqCount];\r
- WEIGHT *NewWeights = new WEIGHT[uSeqCount];\r
-\r
- for (unsigned uSeqIndex = 0; uSeqIndex < m_uSeqCount; ++uSeqIndex)\r
- {\r
- NewSeqs[uSeqIndex] = m_szSeqs[uSeqIndex];\r
- NewNames[uSeqIndex] = m_szNames[uSeqIndex];\r
- NewWeights[uSeqIndex] = m_Weights[uSeqIndex];\r
- }\r
-\r
- for (unsigned uSeqIndex = m_uSeqCount; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- {\r
- char *Seq = new char[uColCount];\r
- NewSeqs[uSeqIndex] = Seq;\r
-#if DEBUG\r
- memset(Seq, '?', uColCount);\r
-#endif\r
- }\r
-\r
- delete[] m_szSeqs;\r
- delete[] m_szNames;\r
- delete[] m_Weights;\r
-\r
- m_szSeqs = NewSeqs;\r
- m_szNames = NewNames;\r
- m_Weights = NewWeights;\r
-\r
- m_uCacheSeqCount = uSeqCount;\r
- m_uCacheSeqLength = uColCount;\r
- m_uColCount = uColCount;\r
- }\r
-\r
-void MSA::FixAlpha()\r
- {\r
- ClearInvalidLetterWarning();\r
- for (unsigned uSeqIndex = 0; uSeqIndex < m_uSeqCount; ++uSeqIndex)\r
- {\r
- for (unsigned uColIndex = 0; uColIndex < m_uColCount; ++uColIndex)\r
- {\r
- char c = GetChar(uSeqIndex, uColIndex);\r
- if (!IsResidueChar(c) && !IsGapChar(c))\r
- {\r
- char w = GetWildcardChar();\r
- // Warning("Invalid letter '%c', replaced by '%c'", c, w);\r
- InvalidLetterWarning(c, w);\r
- SetChar(uSeqIndex, uColIndex, w);\r
- }\r
- }\r
- }\r
- ReportInvalidLetters();\r
- }\r
-\r
-ALPHA MSA::GuessAlpha() const\r
- {\r
-// If at least MIN_NUCLEO_PCT of the first CHAR_COUNT non-gap\r
-// letters belong to the nucleotide alphabet, guess nucleo.\r
-// Otherwise amino.\r
- const unsigned CHAR_COUNT = 100;\r
- const unsigned MIN_NUCLEO_PCT = 95;\r
-\r
- const unsigned uSeqCount = GetSeqCount();\r
- const unsigned uColCount = GetColCount();\r
- if (0 == uSeqCount)\r
- return ALPHA_Amino;\r
-\r
- unsigned uDNACount = 0;\r
- unsigned uRNACount = 0;\r
- unsigned uTotal = 0;\r
- unsigned i = 0;\r
- for (;;)\r
- {\r
- unsigned uSeqIndex = i/uColCount;\r
- if (uSeqIndex >= uSeqCount)\r
- break;\r
- unsigned uColIndex = i%uColCount;\r
- ++i;\r
- char c = GetChar(uSeqIndex, uColIndex);\r
- if (IsGapChar(c))\r
- continue;\r
- if (IsDNA(c))\r
- ++uDNACount;\r
- if (IsRNA(c))\r
- ++uRNACount;\r
- ++uTotal;\r
- if (uTotal >= CHAR_COUNT)\r
- break;\r
- }\r
- if (uTotal != 0 && ((uRNACount*100)/uTotal) >= MIN_NUCLEO_PCT)\r
- return ALPHA_RNA;\r
- if (uTotal != 0 && ((uDNACount*100)/uTotal) >= MIN_NUCLEO_PCT)\r
- return ALPHA_DNA;\r
- return ALPHA_Amino;\r
- }\r
git://source.jalview.org / jabaws.git / blobdiff