#include "muscle.h" #include "textfile.h" #include "msa.h" #include "tree.h" #include "profile.h" #include "objscore.h" bool TreeNeededForWeighting(SEQWEIGHT s) { switch (s) { case SEQWEIGHT_ClustalW: case SEQWEIGHT_ThreeWay: return true; default: return false; } } static ProfPos *ProfileFromMSALocal(MSA &msa, Tree &tree) { const unsigned uSeqCount = msa.GetSeqCount(); for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex) msa.SetSeqId(uSeqIndex, uSeqIndex); if (TreeNeededForWeighting(g_SeqWeight2)) { TreeFromMSA(msa, tree, g_Cluster2, g_Distance2, g_Root1); SetMuscleTree(tree); } return ProfileFromMSA(msa); } void ProfileProfile(MSA &msa1, MSA &msa2, MSA &msaOut) { //ALPHA Alpha = ALPHA_Undefined; //switch (g_SeqType) // { //case SEQTYPE_Auto: // Alpha = msa1.GuessAlpha(); // break; //case SEQTYPE_Protein: // Alpha = ALPHA_Amino; // break; //case SEQTYPE_DNA: // Alpha = ALPHA_DNA; // break; //case SEQTYPE_RNA: // Alpha = ALPHA_RNA; // break; //default: // Quit("Invalid SeqType"); // } //SetAlpha(Alpha); //msa1.FixAlpha(); //msa2.FixAlpha(); unsigned uLength1; unsigned uLength2; uLength1 = msa1.GetColCount(); uLength2 = msa2.GetColCount(); Tree tree1; Tree tree2; ProfPos *Prof1 = ProfileFromMSALocal(msa1, tree1); ProfPos *Prof2 = ProfileFromMSALocal(msa2, tree2); PWPath Path; ProfPos *ProfOut; unsigned uLengthOut; Progress("Aligning profiles"); AlignTwoProfs(Prof1, uLength1, 1.0, Prof2, uLength2, 1.0, Path, &ProfOut, &uLengthOut); Progress("Building output"); AlignTwoMSAsGivenPath(Path, msa1, msa2, msaOut); } // Do profile-profile alignment void Profile() { if (0 == g_pstrFileName1 || 0 == g_pstrFileName2) Quit("-profile needs -in1 and -in2"); SetSeqWeightMethod(g_SeqWeight1); TextFile file1(g_pstrFileName1); TextFile file2(g_pstrFileName2); MSA msa1; MSA msa2; MSA msaOut; Progress("Reading %s", g_pstrFileName1); msa1.FromFile(file1); Progress("%u seqs %u cols", msa1.GetSeqCount(), msa1.GetColCount()); Progress("Reading %s", g_pstrFileName2); msa2.FromFile(file2); Progress("%u seqs %u cols", msa2.GetSeqCount(), msa2.GetColCount()); ALPHA Alpha = ALPHA_Undefined; switch (g_SeqType) { case SEQTYPE_Auto: Alpha = msa1.GuessAlpha(); break; case SEQTYPE_Protein: Alpha = ALPHA_Amino; break; case SEQTYPE_DNA: Alpha = ALPHA_DNA; break; case SEQTYPE_RNA: Alpha = ALPHA_RNA; break; default: Quit("Invalid seq type"); } SetAlpha(Alpha); msa1.FixAlpha(); msa2.FixAlpha(); SetPPScore(); if (ALPHA_DNA == Alpha || ALPHA_RNA == Alpha) SetPPScore(PPSCORE_SPN); const unsigned uSeqCount1 = msa1.GetSeqCount(); const unsigned uSeqCount2 = msa2.GetSeqCount(); const unsigned uSumSeqCount = uSeqCount1 + uSeqCount2; MSA::SetIdCount(uSumSeqCount); ProfileProfile(msa1, msa2, msaOut); Progress("Writing output"); MuscleOutput(msaOut); }