+++ /dev/null
-#include "muscle.h"\r
-#include "msa.h"\r
-#include "seqvect.h"\r
-#include "textfile.h"\r
-\r
-#define MEMDEBUG 0\r
-\r
-#if MEMDEBUG\r
-#include <crtdbg.h>\r
-#endif\r
-\r
-void MUSCLE(SeqVect &v, MSA &msaOut);\r
-\r
-// Append msa2 at the end of msa1\r
-void AppendMSA(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
-\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;\r
- bool bFound = msa2.GetSeqIndex(uId, &uSeqIndex2);\r
- if (bFound)\r
- {\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
- else\r
- {\r
- for (unsigned uColIndex = 0; uColIndex < uColCount2; ++uColIndex)\r
- msa1.SetChar(uSeqIndex, uColCount1 + uColIndex, '-');\r
- }\r
- }\r
- }\r
-\r
-static void SeqFromMSACols(const MSA &msa, unsigned uSeqIndex, unsigned uColFrom,\r
- unsigned uColTo, Seq &s)\r
- {\r
- s.Clear();\r
- s.SetName(msa.GetSeqName(uSeqIndex));\r
- s.SetId(msa.GetSeqId(uSeqIndex));\r
- for (unsigned uColIndex = uColFrom; uColIndex <= uColTo; ++uColIndex)\r
- {\r
- char c = msa.GetChar(uSeqIndex, uColIndex);\r
- if (!IsGapChar(c))\r
- s.AppendChar(c);\r
- }\r
- }\r
-\r
-static void SeqVectFromMSACols(const MSA &msa, unsigned uColFrom, unsigned uColTo,\r
- 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
- SeqFromMSACols(msa, uSeqIndex, uColFrom, uColTo, s);\r
- v.AppendSeq(s);\r
- }\r
- }\r
-\r
-void RefineW(const MSA &msaIn, MSA &msaOut)\r
- {\r
- const unsigned uSeqCount = msaIn.GetSeqCount();\r
- const unsigned uColCount = msaIn.GetColCount();\r
-\r
-// Reserve same nr seqs, 20% more cols\r
- const unsigned uReserveColCount = (uColCount*120)/100;\r
- msaOut.SetSize(uSeqCount, uReserveColCount);\r
-\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- {\r
- msaOut.SetSeqName(uSeqIndex, msaIn.GetSeqName(uSeqIndex));\r
- msaOut.SetSeqId(uSeqIndex, msaIn.GetSeqId(uSeqIndex));\r
- }\r
-\r
- const unsigned uWindowCount = (uColCount + g_uRefineWindow - 1)/g_uRefineWindow;\r
- if (0 == g_uWindowTo)\r
- g_uWindowTo = uWindowCount - 1;\r
-\r
-#if MEMDEBUG\r
- _CrtSetBreakAlloc(1560);\r
-#endif\r
-\r
- if (g_uWindowOffset > 0)\r
- {\r
- MSA msaTmp;\r
- MSAFromColRange(msaIn, 0, g_uWindowOffset, msaOut);\r
- }\r
-\r
- fprintf(stderr, "\n");\r
- for (unsigned uWindowIndex = g_uWindowFrom; uWindowIndex <= g_uWindowTo; ++uWindowIndex)\r
- {\r
- fprintf(stderr, "Window %d of %d \r", uWindowIndex, uWindowCount);\r
- const unsigned uColFrom = g_uWindowOffset + uWindowIndex*g_uRefineWindow;\r
- unsigned uColTo = uColFrom + g_uRefineWindow - 1;\r
- if (uColTo >= uColCount)\r
- uColTo = uColCount - 1;\r
- assert(uColTo >= uColFrom);\r
-\r
- SeqVect v;\r
- SeqVectFromMSACols(msaIn, uColFrom, uColTo, v);\r
-\r
-#if MEMDEBUG\r
- _CrtMemState s1;\r
- _CrtMemCheckpoint(&s1);\r
-#endif\r
-\r
- MSA msaTmp;\r
- MUSCLE(v, msaTmp);\r
- AppendMSA(msaOut, msaTmp);\r
- if (uWindowIndex == g_uSaveWindow)\r
- {\r
- MSA msaInTmp;\r
- unsigned uOutCols = msaOut.GetColCount();\r
- unsigned un = uColTo - uColFrom + 1;\r
- MSAFromColRange(msaIn, uColFrom, un, msaInTmp);\r
-\r
- char fn[256];\r
- sprintf(fn, "win%d_inaln.tmp", uWindowIndex);\r
- TextFile fIn(fn, true);\r
- msaInTmp.ToFile(fIn);\r
-\r
- sprintf(fn, "win%d_inseqs.tmp", uWindowIndex);\r
- TextFile fv(fn, true);\r
- v.ToFile(fv);\r
-\r
- sprintf(fn, "win%d_outaln.tmp", uWindowIndex);\r
- TextFile fOut(fn, true);\r
- msaTmp.ToFile(fOut);\r
- }\r
-\r
-#if MEMDEBUG\r
- void FreeDPMemSPN();\r
- FreeDPMemSPN();\r
-\r
- _CrtMemState s2;\r
- _CrtMemCheckpoint(&s2);\r
-\r
- _CrtMemState s;\r
- _CrtMemDifference(&s, &s1, &s2);\r
-\r
- _CrtMemDumpStatistics(&s);\r
- _CrtMemDumpAllObjectsSince(&s1);\r
- exit(1);\r
-#endif\r
-//#if DEBUG\r
-// AssertMSAEqIgnoreCaseAndGaps(msaInTmp, msaTmp);\r
-//#endif\r
- }\r
- fprintf(stderr, "\n");\r
-\r
-// AssertMSAEqIgnoreCaseAndGaps(msaIn, msaOut);//@@uncomment!\r
- }\r
-\r
-void DoRefineW()\r
- {\r
- SetOutputFileName(g_pstrOutFileName);\r
- SetInputFileName(g_pstrInFileName);\r
- SetStartTime();\r
-\r
- SetMaxIters(g_uMaxIters);\r
- SetSeqWeightMethod(g_SeqWeight1);\r
-\r
- TextFile fileIn(g_pstrInFileName);\r
- MSA msa;\r
- msa.FromFile(fileIn);\r
-\r
- const unsigned uSeqCount = msa.GetSeqCount();\r
- if (0 == uSeqCount)\r
- Quit("No sequences in input file");\r
-\r
- MSA::SetIdCount(uSeqCount);\r
-\r
-// Initialize sequence ids.\r
-// From this point on, ids must somehow propogate from here.\r
- for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
- msa.SetSeqId(uSeqIndex, uSeqIndex);\r
- SetMuscleInputMSA(msa);\r
-\r
- ALPHA Alpha = ALPHA_Undefined;\r
- switch (g_SeqType)\r
- {\r
- case SEQTYPE_Auto:\r
- Alpha = msa.GuessAlpha();\r
- break;\r
-\r
- case SEQTYPE_Protein:\r
- Alpha = ALPHA_Amino;\r
- break;\r
-\r
- case SEQTYPE_DNA:\r
- Alpha = ALPHA_DNA;\r
- break;\r
-\r
- case SEQTYPE_RNA:\r
- Alpha = ALPHA_RNA;\r
- break;\r
-\r
- default:\r
- Quit("Invalid SeqType");\r
- }\r
- SetAlpha(Alpha);\r
- msa.FixAlpha();\r
-\r
- if (ALPHA_DNA == Alpha || ALPHA_RNA == Alpha)\r
- SetPPScore(PPSCORE_SPN);\r
-\r
- MSA msaOut;\r
- RefineW(msa, msaOut);\r
-\r
-// ValidateMuscleIds(msa);\r
-\r
-// TextFile fileOut(g_pstrOutFileName, true);\r
-// msaOut.ToFile(fileOut);\r
- MuscleOutput(msaOut);\r
- }\r
git://source.jalview.org / jabaws.git / blobdiff