+++ /dev/null
-/**
- * Author: Mark Larkin
- *
- * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
- */
-/**
- * Changes:
- *
- * 10-02-07,Nigel Brown(EMBL): changed ifstream to InFileStream to handle
- * cross-platform end-of-lines (for dendrograms, not for help file). Remerged
- * this change 13-04-07.
- *
- * Mark 10-5-2007: Bug fix # 42. call getWeightsForQtLowScore function in
- * QTcalcWeightsForLowScoreSeg instead of getWeightsFromDistMat.
- *
- * Mark 22-5-07, changed the distmatrix to be the size of alignObject.numSeqs
- * Mark Change 20-6-07, added call to calculateMaxLengths()
- * Mark, 3-7-07, Changed getHelp.
- */
-#ifdef HAVE_CONFIG_H
- #include "config.h"
-#endif
-
-#include <iostream>
-#include <fstream>
-#include <cstdio>
-#include "Clustal.h"
-#include "pairwise/FullPairwiseAlign.h"
-#include "pairwise/FastPairwiseAlign.h"
-#include "multipleAlign/MSA.h"
-#include "multipleAlign/LowScoreSegProfile.h"
-#include "multipleAlign/Iteration.h"
-#include "tree/TreeInterface.h"
-#include "general/debuglogObject.h"
-#include "general/statsObject.h"
-#include "alignment/ObjectiveScore.h"
-#include "general/ClustalWResources.h"
-#include "Help.h"
-#include <ctime>
-
-using namespace std;
-
-namespace clustalw
-{
-
-Clustal::Clustal()
-{
-#ifdef WINDOWS
- helpFileName = string("clustalw.hlp");
-#else
- helpFileName = string("clustalw_help");
-#endif
-
- checkTree = true;
- newSeq = 0;
- sequencesMsg = "\nUse the existing GUIDE TREE file, ";
- profile1Msg = "\nUse the existing GUIDE TREE file for Profile 1, " ;
- profile2Msg = "\nUse the existing GUIDE TREE file for Profile 2, ";
-
- newProfile1TreePrompt = "\nEnter name for new GUIDE TREE file for profile 1 [";
- newProfile2TreePrompt = "\nEnter name for new GUIDE TREE file for profile 2 [";
-
- initInterface();
-}
-
-/** ***********************************************************************
- * The function align is used to do a full multiple sequence alignment. *
- **************************************************************************/
-// FIXME: merge doAlignUseOldTree in here
-void Clustal::align(string* phylipName, bool createOutput)
-{
- //time_t start, end;
- //double dif;
- //start = time (NULL);
- //ObjectiveScore score;
- //double _score = score.getSagaScore(&alignmentObj);
- //cout << "SAGA score " << _score << "\n";
-
- string path;
- int count;
- AlignmentOutput alignOutput;
-
- if(userParameters->getEmpty() && userParameters->getMenuFlag())
- {
- utilityObject->error("No sequences in memory. Load sequences first.");
- return;
- }
-
- userParameters->setStructPenalties1(NONE);
- userParameters->setStructPenalties2(NONE);
-
- alignmentObj.clearSecStruct1();
- alignmentObj.clearSecStruct2();
-
- utilityObject->getPath(userParameters->getSeqName(), &path);
-
- if(createOutput && (userParameters->getMenuFlag() || !userParameters->getInteractive()))
- {
- if(!alignOutput.openAlignmentOutput(path))
- {
- return;
- }
- }
- else if(createOutput)
- {
- // We are using clustalQT.
- // Open all the files.
- if(!alignOutput.QTOpenFilesForOutput(QTFileNames))
- {
- return; // could not open the files.
- }
- }
-
- if (userParameters->getSaveParameters())
- {
- userParameters->createParameterOutput();
- }
-
- if(userParameters->getResetAlignmentsNew() || userParameters->getResetAlignmentsAll())
- {
- alignmentObj.resetAlign();
- }
- if(userParameters->getDisplayInfo())
- {
- cout << "Start of Pairwise alignments\n";
- cout << "Aligning...\n";
- }
- if(userParameters->getDNAFlag())
- {
- userParameters->setDNAParams();
- }
- else
- {
- userParameters->setProtParams();
- }
-
- if (statsObject->isEnabled()) {
- statsObject->logInputSeqStats(&alignmentObj);
- }
-
-
- /// STEP 1: PAIRWISE ALIGNMENT TO GENERATE DISTANCE MATRIX
-
- SymMatrix distMat;
- int _numSeqs = alignmentObj.getNumSeqs();
- distMat.ResizeRect(_numSeqs + 1);
-
- PairwiseAlignBase* pairwiseDist;
- if (userParameters->getQuickPairAlign())
- {
- pairwiseDist = new FastPairwiseAlign();
- }
- else
- {
- pairwiseDist = new FullPairwiseAlign();
- }
- // Generate distance matrix!
- pairwiseDist->pairwiseAlign(&alignmentObj, &distMat, 0, _numSeqs, 0, _numSeqs);
- delete pairwiseDist;
-
- bool success = false;
- auto_ptr<AlignmentSteps> progSteps;
- vector<int> seqWeight(_numSeqs + 1);
-
- #if DEBUGFULL
- if(logObject && DEBUGLOG)
- {
- logObject->logMsg("Calling getWeightsAndStepsFromDistMat\n");
- }
- #endif
-
- TreeInterface calcSteps;
- progSteps = calcSteps.getWeightsAndStepsFromDistMat(&seqWeight, &distMat, &alignmentObj,
- 1, _numSeqs, phylipName, &success);
- //cout << "weights and steps calculated!\n";
- //end = time (NULL);
- //dif = difftime(end, start);
- //cout << "It took " << dif << " seconds so Far\n";
-
- if(!success)
- {
- #if DEBUGFULL
- if(logObject && DEBUGLOG)
- {
- logObject->logMsg("Unsuccessful!!!\n");
- }
- #endif
- return;
- }
- #if DEBUGFULL
- if(logObject && DEBUGLOG)
- {
- logObject->logMsg("doing multiSeqAlign\n");
- }
- #endif
- MSA* msaObj = new MSA();
- count = msaObj->multiSeqAlign(&alignmentObj, &distMat, &seqWeight, progSteps.get(), 0);
- delete msaObj;
- //cout << "alignment finished!\n";
-
- //end = time (NULL);
- //dif = difftime(end, start);
- //cout << "It took " << dif << " seconds so Far\n";
-
- if (count <= 0)
- {
- return;
- }
-
- if (userParameters->getMenuFlag())
- {
- cout << "\n\n\n";
- }
-
- /// Do iteration to improve alignment!!!
- if(userParameters->getDoRemoveFirstIteration() == ALIGNMENT)
- {
- //userParameters->setNumIterations(_numSeqs * 2);
- Iteration iterateObj;
- iterateObj.removeFirstIterate(&alignmentObj);
- alignmentObj.calculateMaxLengths(); // Mark Change 20-6-07
- if(userParameters->getDisplayInfo())
- cout << "Finished iteration\n";
- }
-
- if (statsObject->isEnabled()) {
- statsObject->logAlignedSeqStats(&alignmentObj);
- }
-
-
-
- /// STEP 4: OUTPUT THE ALIGNMENT
- if(createOutput)
- {
- alignOutput.createAlignmentOutput(&alignmentObj, 1, _numSeqs);
- }
- (*phylipName) = "";
-
- //end = time (NULL);
- //dif = difftime(end, start);
- //cout << "It took " << dif << " seconds\n";
-}
-
-/** ****************************************************************************
- * The function sequencesAlignToProfile is used to align a set of sequences to *
- * a profile *
- *******************************************************************************/
-void Clustal::sequencesAlignToProfile(string* phylipName)
-{
- string path;
- string treeName;
- bool useTree;
- int i, j, count;
- float dscore;
- bool saveSS2;
- AlignmentOutput alignOutput;
-
- if(userParameters->getProfile1Empty() && userParameters->getMenuFlag())
- {
- utilityObject->error("No profile in memory. Input 1st profile first.\n");
- return;
- }
-
- if(userParameters->getProfile2Empty() && userParameters->getMenuFlag())
- {
- utilityObject->error("No sequences in memory. Input sequences first.\n");
- return;
- }
-
- utilityObject->getPath(userParameters->getProfile2Name(), &path);
-
- if(userParameters->getMenuFlag() || !userParameters->getInteractive())
- {
- if(!alignOutput.openAlignmentOutput(path))
- {
- return;
- }
- }
- else
- {
- // We are using clustalQT.
- // Open all the files.
- if(!alignOutput.QTOpenFilesForOutput(QTFileNames))
- {
- return; // could not open the files.
- }
- }
-
- newSeq = alignmentObj.getProfile1NumSeqs() + 1;
-
- // check for secondary structure information for list of sequences
-
- saveSS2 = userParameters->getUseSS2();
-
- if (userParameters->getStructPenalties2() != NONE && userParameters->getUseSS2() == true
- && (alignmentObj.getNumSeqs() - alignmentObj.getProfile1NumSeqs() > 1))
- {
- if (userParameters->getStructPenalties2() == SECST)
- {
- utilityObject->warning("\n\nWARNING: ignoring secondary structure for a list of sequences\n\n");
- }
- else if (userParameters->getStructPenalties2() == GMASK)
- {
- utilityObject->warning("\n\nWARNING: ignoring gap penalty mask for a list of sequences\n\n");
- }
- userParameters->setUseSS2(false);
- }
-
- DistMatrix distMat;
- int _numSeqs = alignmentObj.getNumSeqs();
- distMat.ResizeRect(_numSeqs + 1);
-
- //
- // Convert to similarities!!!!!!!!
- // This is calcualting the similarities of the sequences in the profile part
- //
- for (i = 1; i <= newSeq; i++)
- {
- for (j = i + 1; j <= newSeq; j++)
- {
- dscore = alignmentObj.countid(i,j);
- distMat(i, j) = ((double)100.0 - (double)dscore)/(double)100.0;
- distMat(j, i) = distMat(i, j);
- }
- }
- //distMat.printArray();
-
- InFileStream _treeFile; //nige
-
- useTree = false;
- //
- // Note put this into a separate function!!!!!
- //
- if (_numSeqs >= 2)
- {
- useTree = useExistingGuideTree(Sequences, phylipName, path);
- }
-
- if (userParameters->getSaveParameters())
- {
- userParameters->createParameterOutput();
- }
-
- if(userParameters->getResetAlignmentsNew() || userParameters->getResetAlignmentsAll())
- {
- alignmentObj.resetProfile2();
- }
- else
- {
- alignmentObj.fixGaps();
- }
-
- int _length = 0;
- if (userParameters->getStructPenalties1() == SECST)
- {
- _length = alignmentObj.getSeqLength(1);
- calcGapPenaltyMask(_length, alignmentObj.getSecStructMask1(),
- alignmentObj.getGapPenaltyMask1());
- }
- if (userParameters->getStructPenalties2() == SECST)
- {
- _length = alignmentObj.getSeqLength(alignmentObj.getProfile1NumSeqs() + 1);
- calcGapPenaltyMask(_length, alignmentObj.getSecStructMask2(),
- alignmentObj.getGapPenaltyMask2());
- }
-
- // PROGRESSIVE ALIGNMENT ALGORITHM //
- vector<int> seqWeights(_numSeqs + 1);
-
- bool success = false;
-
- if (useTree == false) // create the new tree file, if necessary
- {
- if(userParameters->getDisplayInfo())
- {
- cout << "Start of Pairwise alignments\n";
- cout << "Aligning...\n";
- }
-
- if(userParameters->getDNAFlag())
- {
- userParameters->setDNAParams();
- }
- else
- {
- userParameters->setProtParams();
- }
-
- // STEP 1: CALCULATE DISTANCE MATRIX USING PAIRWISE ALIGNMENT //
- PairwiseAlignBase* pairwiseDist;
- if (userParameters->getQuickPairAlign())
- {
- pairwiseDist = new FastPairwiseAlign();
- }
- else
- {
- pairwiseDist = new FullPairwiseAlign();
- }
-
- // Generate distance matrix!
- pairwiseDist->pairwiseAlign(&alignmentObj, &distMat, 0, _numSeqs, newSeq-2, _numSeqs);
- delete pairwiseDist;
-
- if(userParameters->getDisplayInfo())
- cout << "\n\n";
-
- TreeInterface calcSeqWeights;
- calcSeqWeights.getWeightsFromDistMat(&seqWeights, &distMat, &alignmentObj, 1,
- _numSeqs, phylipName, &success);
- }
- else
- {
- TreeInterface calcSeqWeights;
- calcSeqWeights.getWeightsFromGuideTree(&alignmentObj, &distMat, phylipName,
- &seqWeights, 1, _numSeqs, &success);
- }
-
- if(!success)
- {
- return;
- }
- // If it is true, call the function here to get the seqweights etc from it.
-
- // if users request only the guide tree, return
- if (userParameters->getNewTreeFile())
- {
- return;
- }
-
-
- MSA* msaObj = new MSA();
- count = msaObj->seqsAlignToProfile(&alignmentObj, &distMat, &seqWeights, newSeq - 2,
- *phylipName);
- delete msaObj;
-
- userParameters->setUseSS2(saveSS2);
-
- if (count <= 0)
- {
- return;
- }
-
- if (userParameters->getMenuFlag())
- {
- cout << "\n\n\n";
- }
-
- /// STEP 4: OUTPUT THE ALIGNMENT //
- alignOutput.createAlignmentOutput(&alignmentObj, 1, _numSeqs);
-
- (*phylipName) = "";
-}
-
-/** ****************************************************************************
- * The function profileAlign is used to align two profiles *
- *******************************************************************************/
-void Clustal::profileAlign(string* p1TreeName, string* p2TreeName)
-{
- string path;
- //string treeName;
- bool useTree1, useTree2;
- int count, i, j, dscore;
- int _profile1NumSeqs = alignmentObj.getProfile1NumSeqs();
- AlignmentOutput alignOutput;
-
- if(userParameters->getProfile1Empty() || userParameters->getProfile2Empty())
- {
- utilityObject->error("No sequences in memory. Load sequences first.\n\n");
- return;
- }
-
- utilityObject->getPath(userParameters->getProfile1Name(), &path);
-
- if(userParameters->getMenuFlag() || !userParameters->getInteractive())
- {
- if(!alignOutput.openAlignmentOutput(path))
- {
- return;
- }
- }
- else
- {
- // We are using clustalQT.
- // Open all the files.
- if(!alignOutput.QTOpenFilesForOutput(QTFileNames))
- {
- return; // could not open the files.
- }
- }
-
- if(userParameters->getResetAlignmentsNew() || userParameters->getResetAlignmentsAll())
- {
- alignmentObj.resetProfile1();
- alignmentObj.resetProfile2();
- }
- else
- {
- alignmentObj.fixGaps();
- }
-
- // Check if there exists a tree for profile1
- useTree1 = false;
-
- if (_profile1NumSeqs >= 2)
- {
- useTree1 = useExistingGuideTree(Profile1, p1TreeName, path);
- }
-
- // Check if there exists a tree for profile2
- useTree2 = false;
-
- utilityObject->getPath(userParameters->getProfile2Name(), &path);
-
- if (alignmentObj.getNumSeqs() - _profile1NumSeqs >= 2)
- {
- useTree2 = useExistingGuideTree(Profile2, p2TreeName, path);
- }
-
- if (userParameters->getSaveParameters())
- {
- userParameters->createParameterOutput();
- }
- int _length = 0;
- if (userParameters->getStructPenalties1() == SECST)
- {
- _length = alignmentObj.getSeqLength(1);
- calcGapPenaltyMask(_length, alignmentObj.getSecStructMask1(),
- alignmentObj.getGapPenaltyMask1());
- }
-
- if (userParameters->getStructPenalties2() == SECST)
- {
- _length = alignmentObj.getSeqLength(_profile1NumSeqs + 1);
- calcGapPenaltyMask(_length, alignmentObj.getSecStructMask2(),
- alignmentObj.getGapPenaltyMask2());
- }
-
- // Declare the distance matrix
- DistMatrix distMat;
- int _numSeqs = alignmentObj.getNumSeqs();
- distMat.ResizeRect(_numSeqs + 1);
-
- if (useTree1 == false)
- {
- if (_profile1NumSeqs >= 2)
- {
- for (i = 1; i <= _profile1NumSeqs; i++)
- {
- for (j = i + 1; j <= _profile1NumSeqs; j++)
- {
- dscore = static_cast<int>(alignmentObj.countid(i,j));
- distMat(i, j) = (100.0 - dscore)/100.0;
- distMat(j, i) = distMat(i, j);
- }
- }
- utilityObject->getPath(userParameters->getProfile1Name(), &path);
-
- // We need to get the name of the file first because the message is different!
- if(userParameters->getMenuFlag())
- {
- promptForNewGuideTreeName(Profile1, p1TreeName, path);
- }
- else
- {
- string treeName;
- treeName = path + "dnd";
- p1TreeName = new string(treeName);
- }
- }
-
- if (useTree2 == false)
- {
- if(_numSeqs - _profile1NumSeqs >= 2)
- {
- for (i = 1 + _profile1NumSeqs; i <= _numSeqs; i++)
- {
- for (j = i + 1; j <= _numSeqs; j++)
- {
- dscore = static_cast<int>(alignmentObj.countid(i,j));
- distMat(i, j) = (100.0 - dscore) / 100.0;
- distMat(j, i) = distMat(i, j);
- }
- }
-
- utilityObject->getPath(userParameters->getProfile2Name(), &path);
-
- if(userParameters->getMenuFlag())
- {
- promptForNewGuideTreeName(Profile2, p2TreeName, path);
- }
- else
- {
- string treeName;
- treeName = path + "dnd";
- p2TreeName = new string(treeName);
- }
-
- }
- }
- }
-
- bool success = false;
- vector<int> prof1Weight, prof2Weight;
- prof1Weight.resize(_profile1NumSeqs);
- prof2Weight.resize(_numSeqs);
- TreeInterface tree;
- tree.getWeightsForProfileAlign(&alignmentObj, &distMat, p1TreeName, &prof1Weight,
- p2TreeName, &prof2Weight, _numSeqs, _profile1NumSeqs,
- useTree1, useTree2, &success);
- if(!success)
- {
- return;
- }
-
- // do an initial alignment to get the pairwise identities between the two
- // profiles - used to set parameters for the final alignment
- MSA* msaObj = new MSA();
-
- alignmentObj.resetProfile1();
- alignmentObj.resetProfile2();
-
- count = msaObj->doProfileAlign(&alignmentObj, &distMat, &prof1Weight, &prof2Weight);
- delete msaObj;
-
- if (count == 0)
- {
- return;
- }
- if(userParameters->getMenuFlag())
- {
- cout << "\n\n\n";
- }
-
- alignOutput.createAlignmentOutput(&alignmentObj, 1, _numSeqs);
-
- (*p1TreeName) = "";
- (*p2TreeName) = "";
-}
-
-void Clustal::doGuideTreeOnly(string* phylipName)
-{
- string path;
- if(userParameters->getEmpty())
- {
- utilityObject->error("No sequences in memory. Load sequences first.\n");
- return;
- }
-
- userParameters->setStructPenalties1(NONE);
- userParameters->setStructPenalties2(NONE);
-
- alignmentObj.clearSecStruct1();
- alignmentObj.clearSecStruct2();
-
- if(userParameters->getResetAlignmentsNew() || userParameters->getResetAlignmentsAll())
- {
- alignmentObj.resetAlign();
- }
-
- utilityObject->getPath(userParameters->getSeqName(), &path);
- int _numSeqs = alignmentObj.getNumSeqs();
-
- if (_numSeqs < 2)
- {
- utilityObject->error("Less than 2 sequences in memory. Phylogenetic tree cannot be built.\n");
- return;
- }
-
- if (userParameters->getSaveParameters())
- {
- userParameters->createParameterOutput();
- }
-
- if(userParameters->getDisplayInfo())
- {
- cout << "Start of Pairwise alignments\n";
- cout << "Aligning...\n";
- }
-
- if(userParameters->getDNAFlag())
- {
- userParameters->setDNAParams();
- }
- else
- {
- userParameters->setProtParams();
- }
-
- ///STEP 1: PAIRWISE ALIGNMENT TO GENERATE DISTANCE MATRIX //
-
- DistMatrix distMat;
- distMat.ResizeRect(_numSeqs + 1);
-
- PairwiseAlignBase* pairwiseDist;
- if (userParameters->getQuickPairAlign())
- {
- pairwiseDist = new FastPairwiseAlign();
- }
- else
- {
- pairwiseDist = new FullPairwiseAlign();
- }
- // Generate distance matrix!
- pairwiseDist->pairwiseAlign(&alignmentObj, &distMat, 0, _numSeqs, 0, _numSeqs);
- delete pairwiseDist;
-
- /* AW DEBUG
- fprintf(stdout, "\nDEBUG: distance matrix following...\n");
- distMat.printArray();
- */
-
- bool success = false;
- TreeInterface tree;
- tree.generateTreeFromDistMat(&distMat, &alignmentObj, 1, _numSeqs, phylipName, &success);
-
- if(userParameters->getResetAlignmentsNew() || userParameters->getResetAlignmentsAll())
- {
- alignmentObj.resetAlign();
- }
-
- (*phylipName) = "";
-}
-
-
-// FIXME this is to 90% identical with align(), please merge
-void Clustal::doAlignUseOldTree(string* phylipName)
-{
- string path;
- int count;
- AlignmentOutput alignOutput;
-
- if(userParameters->getEmpty())
- {
- utilityObject->error("No sequences in memory. Load sequences first.\n");
- return;
- }
-
- userParameters->setStructPenalties1(NONE);
- userParameters->setStructPenalties2(NONE);
-
- alignmentObj.clearSecStruct1();
- alignmentObj.clearSecStruct2();
-
- utilityObject->getPath(userParameters->getSeqName(), &path);
-
- if(userParameters->getMenuFlag() || !userParameters->getInteractive())
- {
- if(!alignOutput.openAlignmentOutput(path))
- {
- return;
- }
- }
- else
- {
- // We are using clustalQT.
- // Open all the files.
- if(!alignOutput.QTOpenFilesForOutput(QTFileNames))
- {
- return; // could not open the files.
- }
- }
-
- if(userParameters->getResetAlignmentsNew() || userParameters->getResetAlignmentsAll())
- {
- alignmentObj.resetAlign();
- }
-
- int _numSeqs = alignmentObj.getNumSeqs();
- DistMatrix distMat;
- distMat.ResizeRect(_numSeqs + 1);
- utilityObject->getPath(userParameters->getSeqName(), &path);
-
- if (_numSeqs >= 2)
- {
-
- if(userParameters->getMenuFlag())
- {
- phylipName = new string(path);
- *phylipName = *phylipName + "dnd";
-
- string message, answer;
- message = "\nEnter a name for the guide tree file [" + *phylipName + "]";
- utilityObject->getStr(message, answer);
-
- if(!answer.empty())
- {
- phylipName = new string(answer);
- }
- }
-
- if(userParameters->getMenuFlag() || !userParameters->getInteractive())
- {
- InFileStream _treeFile; //nige
- _treeFile.open(phylipName->c_str());
- if(!_treeFile.is_open())
- {
- utilityObject->error("Cannot open tree file [%s]\n", phylipName->c_str());
- return;
- }
- _treeFile.close();
- }
- }
- else
- {
- if(userParameters->getDisplayInfo())
- {
- cout << "Start of Pairwise alignments\n";
- cout << "Aligning...\n";
- }
- if(userParameters->getDNAFlag())
- {
- userParameters->setDNAParams();
- }
- else
- {
- userParameters->setProtParams();
- }
-
- PairwiseAlignBase* pairwiseDist;
- if (userParameters->getQuickPairAlign())
- {
- pairwiseDist = new FastPairwiseAlign();
- }
- else
- {
- pairwiseDist = new FullPairwiseAlign();
- }
- // Generate distance matrix!
- pairwiseDist->pairwiseAlign(&alignmentObj, &distMat, 0, _numSeqs, 0, _numSeqs);
- delete pairwiseDist;
- }
-
- if (userParameters->getSaveParameters())
- {
- userParameters->createParameterOutput();
- }
- auto_ptr<AlignmentSteps> progSteps;
- vector<int> seqWeights(_numSeqs + 1);
- bool success = false;
- TreeInterface tree;
- progSteps = tree.getWeightsAndStepsFromTree(&alignmentObj, &distMat, phylipName,
- &seqWeights, 1, _numSeqs, &success);
-
- if(!success)
- {
- return;
- }
-
- MSA* msaObj = new MSA();
- count = msaObj->multiSeqAlign(&alignmentObj, &distMat, &seqWeights, progSteps.get(), 0);
- delete msaObj;
-
- if (count <= 0)
- {
- return;
- }
-
- if (userParameters->getMenuFlag())
- {
- cout << "\n\n\n";
- }
-
- // same as in align()
- if(userParameters->getDoRemoveFirstIteration() == ALIGNMENT)
- {
- //userParameters->setNumIterations(_numSeqs * 2);
- Iteration iterateObj;
- iterateObj.removeFirstIterate(&alignmentObj);
- alignmentObj.calculateMaxLengths(); // Mark Change 20-6-07
- if(userParameters->getDisplayInfo())
- cout << "Finished iteration\n";
- }
-
-
- alignOutput.createAlignmentOutput(&alignmentObj, 1, _numSeqs);
-
- phylipName = new string("");
-}
-
-
-
-void Clustal::getFullHelp()
-{
- vector<string> markers;
- Help myhelp;
- bool showtitle = true;
-
- markers = myhelp.ListSectionMarkers();
- for (unsigned int i=0; i<markers.size(); i++) {
- string m = markers[i];
- getHelp(m, showtitle);
- }
-}
-
-void Clustal::getHelp(char helpPointer, bool printTitle)
-{
- string s = "";
- s += helpPointer;
- Clustal::getHelp(s, printTitle);
-}
-
-
-
-/*
- * Andreas Wilm (UCD): edited to support new help system (separate
- * file now which is compiled in)
- *
- * Author?
- * The clustal help file is called clustalw_help. Should be in the same
- * directory. I have changed it to a C++ style implementation. I am taking
- * out support for VMS until later. We will see if we are still supporting that platform.
- */
-void Clustal::getHelp(string helpPointer, bool printTitle)
-{
- Help myhelp;
- string helpString;
-
- helpString = myhelp.GetSection(helpPointer);
-
- if (printTitle) {
- helpString = "\n\n>> HELP " +
- helpPointer +
- " << " +
- myhelp.GetSectionTitle(helpPointer) +
- "\n\n" +
- helpString;
- }
-
-
- if(! userParameters->getMenuFlag())
- {
- cout << helpString;
- }
- else
- {
- string::size_type lastPos = 0;
- string::size_type pos = helpString.find_first_of("\n", lastPos);
- int nlines = 0;
-
- while (pos != string::npos)
- {
- cout << helpString.substr(lastPos, pos - lastPos);
- nlines++;
-
- if(nlines >= PAGE_LEN)
- {
- char tempChar;
- cout << "\nPress [RETURN] to continue or X to stop ";
- cin.get(tempChar);
- if(toupper(tempChar) == 'X')
- {
- return;
- }
- else
- {
- nlines = 0;
- }
- }
-
- lastPos = pos; //helpString.find_first_not_of("\n", pos);
- pos = helpString.find_first_of("\n", lastPos+1);
- //cerr << "DEBUG: pos=" << pos << " lastPos=" << lastPos << "/" << helpString.length() << "\n";
- }
- }
-}
-
-
-
-/**
- * The wrap functions will be used with interface classes to perform the tasks
- * that were previously done there.
- */
-int Clustal::sequenceInput(bool append, string *offendingSeq)
-{
- int code;
- // If we are not appending, we need to clear the Alignment object.
- if(!append)
- {
- alignmentObj.clearAlignment();
- }
-
- FileReader readSeqFile;
- code = readSeqFile.seqInput(&alignmentObj, append, offendingSeq);
- return code;
-}
-
-/**
- * profile1Input is a wrapper function for the profileInput. This is because the
- * other classes dont have access to FileReader
- */
-int Clustal::profile1Input(string profile1Name)
-{
- int code;
- // I need to clear out the Alignment object.
- alignmentObj.clearAlignment();
- userParameters->setProfileNum(1);
-
- userParameters->setSeqName(profile1Name);
- userParameters->setProfile1Name(profile1Name);
-
- FileReader readProfileFile;
- code = readProfileFile.profileInput(&alignmentObj);
-
- // If we are using the commandline check if there are seqs!
- if(!userParameters->getInteractive())
- {
- // AW: FIXME code should be handled higher up the stack and check all codes
- // also, shouldnt we use utilityObject->error()?
- if(code != OK)
- {
- if (code==NOSEQUENCESINFILE)
- cerr << "ERROR: There are no sequences in profile2 file." << std::endl;
- else if (code==ALLNAMESNOTDIFFERENT)
- cerr << "ERROR: Not all sequence names are different" << std::endl;
- else
- cerr << "ERROR: Unhandled error code (" << code << ") returned from profileInput.\n";
- // AW: should we really exit here? What if called from clustalx?
- exit(2);
- }
- }
-
- return code;
-}
-
-/**
- * profile2Input is a wrapper function for the profileInput. This is because the
- * other classes dont have access to FileReader
- */
-int Clustal::profile2Input(string profile2Name)
-{
- int code;
-
- if(userParameters->getProfileNum() == 2)
- {
- // Remove the sequences from the previous one.
- int numSeqsProfile1 = alignmentObj.getProfile1NumSeqs();
- alignmentObj.resizeSeqArray(numSeqsProfile1 + 1);
- // Clear anything from profile2 in alignment.
- alignmentObj.clearSecStruct2();
- }
- userParameters->setProfileNum(2);
-
- userParameters->setSeqName(profile2Name);
- userParameters->setProfile2Name(profile2Name);
-
- FileReader readProfileFile;
- code = readProfileFile.profileInput(&alignmentObj);
-
- if(!userParameters->getInteractive())
- {
- // AW: FIXME code should be handled higher up the stack and check all codes
- // also, shouldnt we use utilityObject->error()?
- if(code != OK) {
- if (code==NOSEQUENCESINFILE)
- cerr << "ERROR: There are no sequences in profile2 file." << std::endl;
- else if (code==ALLNAMESNOTDIFFERENT)
- cerr << "ERROR: Not all sequence names are different" << std::endl;
- else
- cerr << "ERROR: Unhandled error code (" << code << ") returned from profileInput.\n";
- // AW: should we really exit here? What if called from clustalx?
- exit(2);
- }
- }
- return code;
-}
-
-
-/**
- * The function commandline_readseq is called by the command line to
- * read in the sequences. It also prints out the names. This was
- * previously done by the command line parser.
- */
-int Clustal::commandLineReadSeq(int firstSeq)
-{
- // Clear the alignment, although obviously there shouldnt be anything in it.
- alignmentObj.clearAlignment();
- userParameters->setProfileNum(0);
- int code = 0;
- string offendingSeq;
- FileReader readInputFile;
- code = readInputFile.readSeqs(&alignmentObj, firstSeq, &offendingSeq);
-
- if(code != OK)
- {
- if(code == CANNOTOPENFILE)
- {
- utilityObject->error("Cannot open input file. No alignment!\n");
- }
- else if(code == NOSEQUENCESINFILE)
- {
- utilityObject->error("No sequences in file. No alignment!\n");
- }
- else if(code == ALLNAMESNOTDIFFERENT)
- {
- utilityObject->error("Multiple sequences found with same name (found %s at least twice)!", offendingSeq.c_str());
- }
- else if(code == EMPTYSEQUENCE)
- {
- utilityObject->error("Empty sequences found: %s\n", offendingSeq.c_str());
- }
- else if(code == SEQUENCETOOBIG)
- {
- utilityObject->error("Sequence(s) too big: %s\n", offendingSeq.c_str());
- }
- else if(code == BADFORMAT)
- {
- utilityObject->error("Sequences are badly formatted!\n");
- }
- else
- {
- utilityObject->error("\nThere was a problem reading in the file. No alignment!\n");
- }
- exit(-1);
- }
-
- alignmentObj.printSequencesAddedInfo();
-
- userParameters->setEmpty(false);
- return code;
-}
-
-/*
- * The function outputNow is used to output the alignment. It can be called by the
- * menu or command line parser.
- */
-void Clustal::outputNow()
-{
- if(alignmentObj.getNumSeqs() > 0)
- {
- string path = "";
- if(!userParameters->getMenuFlag())
- {
- string _seqName = userParameters->getSeqName();
- utilityObject->getPath(_seqName, &path);
- }
- AlignmentOutput alignmentOutput;
- alignmentOutput.openAlignmentOutput(path);
-
- alignmentOutput.createAlignmentOutput(&alignmentObj, 1, alignmentObj.getNumSeqs());
- }
- else
- {
- utilityObject->error("No sequences have been loaded\n");
- }
-}
-
-void Clustal::phylogeneticTree(string* phylipName, string* clustalName, string* distName,
- string* nexusName, string pimName)
-{
- TreeNames treeNames;
- treeNames.clustalName = *clustalName;
- treeNames.distName = *distName;
- treeNames.nexusName = *nexusName;
- treeNames.phylipName = *phylipName;
- treeNames.pimName = pimName;
- TreeInterface tree;
- tree.treeFromAlignment(&treeNames, &alignmentObj);
-}
-
-void Clustal::bootstrapTree(string* phylipName, string* clustalName, string* nexusName)
-{
- TreeNames treeNames;
- treeNames.clustalName = *clustalName;
- treeNames.nexusName = *nexusName;
- treeNames.phylipName = *phylipName;
- TreeInterface tree;
- tree.bootstrapTree(&treeNames, &alignmentObj);
-}
-
-void Clustal::initInterface()
-{
- userParameters->setEmpty(true);
- userParameters->setProfile1Empty(true);
- userParameters->setProfile2Empty(true);
-}
-
-/**
- * The function calcGapPenaltyMask is used to calculate the gapMask from the secondary
- * structure information.
- */
-void Clustal::calcGapPenaltyMask(int prfLength, vector<char>* mask, vector<char>* gapMask)
-{
- int i,j;
-
- vector<char> structMask;
- structMask.resize(prfLength + 1);
-
- int _helixEndPlus = userParameters->getHelixEndPlus();
- int _helixEndMinus = userParameters->getHelixEndMinus();
- int _strandEndPlus = userParameters->getStrandEndPlus();
- int _strandEndMinus = userParameters->getStrandEndMinus();
-
- i = 0;
- while (i < prfLength)
- {
- if (tolower((*mask)[i]) == 'a' || (*mask)[i] == '$')
- {
- for (j = -_helixEndPlus; j < 0; j++)
- {
- if(i + j < (int)structMask.size())
- {
- if ((i + j >= 0) && (tolower(structMask[i + j]) != 'a') &&
- (tolower(structMask[i + j]) != 'b'))
- {
- structMask[i + j] = 'a';
- }
- }
- }
- for (j = 0; j < _helixEndMinus; j++)
- {
- if (i + j >= prfLength || (tolower((*mask)[i + j]) != 'a'
- && (*mask)[i + j] != '$'))
- {
- break;
- }
- structMask[i + j] = 'a';
- }
- i += j;
- while (tolower((*mask)[i]) == 'a' || (*mask)[i] == '$')
- {
- if (i >= prfLength)
- {
- break;
- }
- if ((*mask)[i] == '$')
- {
- structMask[i] = 'A';
- i++;
- break;
- }
- else
- {
- structMask[i] = (*mask)[i];
- }
- i++;
- }
- for (j = 0; j < _helixEndMinus; j++)
- {
- if ((i - j - 1 >= 0) && (tolower((*mask)[i - j - 1]) == 'a'
- || (*mask)[i - j - 1] == '$'))
- {
- structMask[i - j - 1] = 'a';
- }
- }
- for (j = 0; j < _helixEndPlus; j++)
- {
- if (i + j >= prfLength)
- {
- break;
- }
- structMask[i+j] = 'a';
- }
- }
- else if (tolower((*mask)[i]) == 'b' || (*mask)[i] == '%')
- {
- for (j = -_strandEndPlus; j < 0; j++)
- {
- if ((i + j >= 0) && (tolower(structMask[i + j]) != 'a')
- && (tolower(structMask[i + j]) != 'b'))
- {
- structMask[i + j] = 'b';
- }
- }
- for (j = 0; j < _strandEndPlus; j++)
- {
- if (i + j >= prfLength || (tolower((*mask)[i + j]) != 'b'
- && (*mask)[i + j] != '%'))
- {
- break;
- }
- structMask[i+j] = 'b';
- }
- i += j;
- while (tolower((*mask)[i]) == 'b' || (*mask)[i] == '%')
- {
- if (i >= prfLength)
- {
- break;
- }
- if ((*mask)[i] == '%')
- {
- structMask[i] = 'B';
- i++;
- break;
- }
- else
- {
- structMask[i] = (*mask)[i];
- }
- i++;
- }
- for (j = 0; j < _strandEndMinus; j++)
- {
- if ((i-j-1>=0) && (tolower((*mask)[i-j-1]) == 'b' || (*mask)[i-j-1] == '%'))
- {
- structMask[i - j - 1] = 'b';
- }
- }
- for (j = 0; j < _strandEndPlus; j++)
- {
- if (i + j >= prfLength)
- {
- break;
- }
- structMask[i+j] = 'b';
- }
- }
- else
- {
- i++;
- }
- }
-
- for(i = 0; i < prfLength;i++)
- {
- switch (structMask[i])
- {
- case 'A':
- (*gapMask)[i] = userParameters->getHelixPenalty() + '0';
- break;
- case 'a':
- (*gapMask)[i] = userParameters->getHelixEndPenalty() + '0';
- break;
- case 'B':
- (*gapMask)[i] = userParameters->getStrandPenalty() +'0';
- break;
- case 'b':
- (*gapMask)[i] = userParameters->getStrandEndPenalty() + '0';
- break;
- default:
- (*gapMask)[i] = userParameters->getLoopPenalty() + '0';
- break;
- }
- }
-}
-
-
-/**
- * The function QTcalcLowScoreSegments is used to calculate the residues in the sequences
- * that score badly.
- * @param firstSeq first seq in the alignment or profile
- * @param nSeqs the number of sequences
- * @param nCols the length of the longest seq
- * @param seqWeight a vector to hold the sequence weights
- * @param lowScoreRes
- * @param seqWeightCalculated
- */
-void Clustal::QTcalcLowScoreSegments(LowScoreSegParams* params)
-{
- int i, j;
- float sum, prevSum;
- float gscale;
- vector<int> weight;
- int sweight;
- vector<int> gaps;
- int matrix[NUMRES][NUMRES];
- vector<float> fsum;
- vector<float> bsum;
- vector<float> pscore;
- int _maxAA = userParameters->getMaxAA();
-
- // STEP 1: Calculate the sequence weights
- QTcalcWeightsForLowScoreSeg(params);
-
- subMatrix->getQTMatrixForLowScoreSeg(matrix);
-
- const SeqArray* seqArray = alignmentObj.getSeqArray();
-
- gaps.resize(params->nCols + 1);
-
- for (j = 1; j <= params->nCols; j++)
- {
- gaps[j - 1] = 0;
- for(i = params->firstSeq + 1; i < params->firstSeq + params->nSeqs; i++)
- {
- if (j < alignmentObj.getSeqLength(i))
- {
- if (((*seqArray)[i][j] < 0) || ((*seqArray)[i][j] > _maxAA))
- {
- gaps[j-1]++;
- }
- }
- }
- }
-
- // STEP 2: Calculate the profile
- LowScoreSegProfile lowScoreProfile(params->nCols, params->firstSeq,
- params->firstSeq + params->nSeqs);
-
- weight.resize(params->firstSeq + params->nSeqs + 1);
-
- for(i = params->firstSeq;i < params->firstSeq + params->nSeqs;i++)
- {
- weight[i] = (*(params->seqWeight))[i - params->firstSeq];
- }
-
- lowScoreProfile.calcLowScoreSegProfile(seqArray, matrix, &weight);
-
- const SeqArray* profile = lowScoreProfile.getProfilePtr();
-
- sweight = 0;
- for(i = params->firstSeq; i < params->firstSeq + params->nSeqs; i++)
- {
- sweight += weight[i];
- }
-
- //Now, use the profile scores to mark segments of each sequence which score badly.
-
- fsum.resize(params->nCols + 2);
- bsum.resize(params->nCols + 2);
- pscore.resize(params->nCols + 2);
-
- for(i = params->firstSeq + 1; i < params->firstSeq + params->nSeqs + 1; i++)
- {
- // In a forward phase, sum the profile scores. Mark negative sums as exceptions.
- //If the sum is positive, then it gets reset to 0.
- sum = 0.0;
- for(j = 1; j <= alignmentObj.getSeqLength(i); j++)
- {
- gscale = (float)(params->nSeqs - gaps[j - 1]) / (float)params->nSeqs;
- if((*seqArray)[i][j] < 0 || (*seqArray)[i][j] >= _maxAA)
- {
- pscore[j - 1] = 0.0;
- sum = 0.0;
- }
- else
- {
- pscore[j-1]=((*profile)[j][(*seqArray)[i][j]]- weight[i - 1] *
- matrix[(*seqArray)[i][j]][(*seqArray)[i][j]]) * gscale / sweight;
- }
- sum += pscore[j - 1];
- if(sum > 0.0)
- {
- sum = 0.0;
- }
- fsum[j - 1] = sum;
- }
-// trim off any positive scoring residues from the end of the segments
- prevSum = 0;
- for(j = alignmentObj.getSeqLength(i) - 1; j >= 0; j--)
- {
- if(prevSum >= 0.0 && fsum[j] < 0.0 && pscore[j] >= 0.0)
- {
- fsum[j] = 0.0;
- }
- prevSum = fsum[j];
- }
-
-// Now, in a backward phase, do the same summing process.
- sum = 0.0;
- for(j = alignmentObj.getSeqLength(i); j >= 1; j--)
- {
- if((*seqArray)[i][j] < 0 || (*seqArray)[i][j] >= _maxAA)
- {
- sum = 0;
- }
- else
- {
- sum += pscore[j - 1];
- }
- if(sum > 0.0)
- {
- sum = 0.0;
- }
- bsum[j - 1] = sum;
- }
-// trim off any positive scoring residues from the start of the segments
- prevSum = 0;
- for(j = 0; j < alignmentObj.getSeqLength(i); j++)
- {
- if(prevSum >= 0.0 && bsum[j] < 0.0 && pscore[j] >= 0.0)
- {
- bsum[j] = 0.0;
- }
- prevSum = bsum[j];
- }
-//Mark residues as exceptions if they score negative in the forward AND backward directions.
- for(j = 1; j <= alignmentObj.getSeqLength(i); j++)
- {
- if(fsum[j - 1] < 0.0 && bsum[j - 1] < 0.0)
- {
- if((*seqArray)[i][j] >= 0 && (*seqArray)[i][j] < _maxAA)
- {
- (*(params->lowScoreRes))[i - params->firstSeq - 1][j - 1] = -1;
- }
- }
- }
- }
-
-// Finally, apply the length cutoff to the segments - removing segments shorter
-//than the cutoff
-
- QTremoveShortSegments(params);
-
-}
-
-void Clustal::QTremoveShortSegments(LowScoreSegParams* params)
-{
- int i,j,k,start;
- //panel_data data;
-
- //GetPanelExtra(p,&data);
- if(params->nSeqs <= 0)
- return;
-
-// Reset all the exceptions - a value of 1 indicates an exception that
-// will be displayed. A value of -1 is used to remember exceptions that
-// are temporarily hidden in the display
-
- for(i = 0; i < params->nSeqs; i++)
- {
- for(j = 0; j < params->nCols; j++)
- {
- if((*(params->lowScoreRes))[i][j] == -1)
- {
- (*(params->lowScoreRes))[i][j] = 1;
- }
- }
- }
-
- for(i = 0; i < params->nSeqs; i++)
- {
- start = -1;
- for(j = 0; j <= params->nCols; j++)
- {
- if(start == -1)
- {
- if((*(params->lowScoreRes))[i][j] == 1)
- start = j;
- }
- else
- {
- if(j == params->nCols || (*(params->lowScoreRes))[i][j] == 0)
- {
- if(j - start < userParameters->getQTminLenLowScoreSegment())
- {
- for(k = start; k < j; k++)
- {
- (*(params->lowScoreRes))[i][k] = -1;
- }
- }
- start = -1;
- }
- }
- }
- }
-}
-
-/**
- * Change: Mark 22-5-07, changed the distmatrix to be the size of alignObject.numSeqs
- */
-void Clustal::QTcalcWeightsForLowScoreSeg(LowScoreSegParams* params)
-{
- int i, j;
- vector<int> weight;
- float dscore;
- DistMatrix distMat(alignmentObj.getNumSeqs() + 1); // Mark: changed size
- // Aw potential trouble here: what if we don't have write
- // permission to current directory?
-#ifdef UNIX
- char treeName[FILENAMELEN]=".score.ph";
-#else
- char treeName[FILENAMELEN]="tmp.ph";
-#endif
-
- if(params->nSeqs <= 0)
- {
- return;
- }
-
-// if sequence weights have been calculated before - don't bother
-//doing it again (it takes too long). data.seqweight is set to NULL when
-// new sequences are loaded. /
- if(params->seqWeightCalculated == true)
- {
- return;
- }
-
- utilityObject->info("Calculating sequence weights...");
-
- // count pairwise percent identities to make a phylogenetic tree //
- if(params->nSeqs >= 2)
- { //i=firstSeq + 1; i <= firstSeq + nSeqs
- for (i = params->firstSeq + 1; i <= params->firstSeq + params->nSeqs; i++)
- { //j=i + 1; i <= firstSeq + nSeqs
- for (j = i + 1; j <= params->firstSeq + params->nSeqs; j++) // Mark 22-5-07
- {
- dscore = alignmentObj.countid(i, j); // Mark 22-5-07
- distMat(i, j) = (100.0 - dscore) / 100.0;
- distMat(j, i) = distMat(i, j);
- }
- }
- string name = string(treeName);
- bool success = false;
- weight.resize(params->firstSeq + params->nSeqs + 1);
- TreeInterface tree;
- tree.getWeightsForQtLowScore(&weight, &distMat, &alignmentObj, params->firstSeq + 1,
- params->nSeqs, &name, &success); // Mark change 10-5-07
- if(!success)
- {
- return;
- }
-
- for(i = params->firstSeq;i < params->firstSeq + params->nSeqs;i++)
- {
- (*(params->seqWeight))[i - params->firstSeq] = weight[i];
- }
-
- utilityObject->info("Done.");
- }
-}
-
-void Clustal::QTSetFileNamesForOutput(AlignmentFileNames fileNames)
-{
- QTFileNames = fileNames;
-}
-
-bool Clustal::QTRealignSelectedRange(AlignmentFileNames fileNames, int beginPos, int endPos, bool realignEndGapPen)
-{
- bool alignEndGapPen = userParameters->getEndGapPenalties();
-
- Alignment saveOldAlign = alignmentObj; // Take a copy of it. Note provided copy
- // constructor is ok.
- bool ok;
- ok = alignmentObj.removeAllOutsideRange(beginPos, endPos);
-
- if(!ok)
- {
- alignmentObj = saveOldAlign;
- return false;
- }
- // Temporarily set the alignment output to be input
- int saveOutOrder = userParameters->getOutputOrder();
- userParameters->setOutputOrder(INPUT);
-
- //set end gap penalties to be realignEndGapPenalties
- userParameters->setEndGapPenalties(realignEndGapPen);
- //do the alignment
- if(alignmentObj.getNumSeqs() <= 0)
- {
- alignmentObj = saveOldAlign;
- return false;
- }
- QTSetFileNamesForOutput(fileNames);
- string phylipName = fileNames.treeFile;
- align(&phylipName, false);
-
- userParameters->setOutputOrder(saveOutOrder);
-
- // reset the end gap penalties
- userParameters->setEndGapPenalties(alignEndGapPen);
-
- // remove postions that only contain gaps
- int nSeqs = alignmentObj.getNumSeqs();
- alignmentObj.removeAllGapOnlyColumns(1, nSeqs, 0);
-
- // save it to a temporary area.
- SeqArray realignedArea = *(alignmentObj.getSeqArray());
-
- // Paste it back into the original alignment.
- alignmentObj = saveOldAlign;
- bool result;
- result = alignmentObj.updateRealignedRange(realignedArea, beginPos, endPos);
- if(result == false)
- {
- utilityObject->error("something went wrong while updating the realigned range\n");
- }
-
- // output the alignments
- AlignmentOutput alignOutput;
- if(!alignOutput.QTOpenFilesForOutput(QTFileNames))
- {
- return false; // could not open the files.
- }
-
- alignOutput.createAlignmentOutput(&alignmentObj, 1, nSeqs);
-
- return true;
-}
-
-void Clustal::test()
-{
- cout << "RUNNING TEST\n";
- AlignmentOutput alignOutput;
- string path;
- utilityObject->getPath(userParameters->getSeqName(), &path);
-
- if(!alignOutput.openAlignmentOutput(path))
- {
- cerr << "could not open the file\n";
- return;
- }
-
- vector<int> selected;
- int nSeqs = alignmentObj.getNumSeqs();
- selected.resize(nSeqs + 1, 0);
- selected[9] = 1;
- selected[10] = 1;
- //selected[1] = 1;
- alignmentObj.removeGapOnlyColsFromSelectedSeqs(&selected);
- alignOutput.createAlignmentOutput(&alignmentObj, 1, nSeqs);
-}
-
-/**
- * This function is used to ask the user if they would like to use an existing guide tree.
- * Note: It expects that phylipName actually points to a string that has been allocated.
- */
-bool Clustal::useExistingGuideTree(int type, string* phylipName, const string& path)
-{
- bool useTree = false;
- string treeName;
- InFileStream _treeFile; //nige
- bool paramUseTree;
-
- string* ptrToMsg;
- if(type == Sequences)
- {
- ptrToMsg = &sequencesMsg;
- paramUseTree = userParameters->getUseTreeFile();
- }
- else if(type == Profile1)
- {
- ptrToMsg = &profile1Msg;
- paramUseTree = userParameters->getUseTree1File();;
- }
- else if(type == Profile2)
- {
- ptrToMsg = &profile2Msg;
- paramUseTree = userParameters->getUseTree2File();;
- }
- else
- {
- ptrToMsg = &sequencesMsg;
- paramUseTree = userParameters->getUseTreeFile();
- }
-
- if (checkTree && userParameters->getMenuFlag())
- {
- treeName = path + "dnd";
-
- _treeFile.open(treeName.c_str());
- _treeFile.seekg(0, std::ios::beg);
-
- if(_treeFile.is_open())
- {
- string message = *ptrToMsg + treeName + " (y/n) ? [y]";
- string answer;
- utilityObject->getStr(message, answer);
-
- if(answer[0] != 'n' && answer[0] != 'N')
- {
- if(!phylipName)
- {
- phylipName = new string(treeName);
- }
- else
- {
- *phylipName = treeName;
- }
- useTree = true;
- }
- _treeFile.close();
- }
- }
- else if (!userParameters->getMenuFlag() && paramUseTree)
- {
- useTree = true;
- }
-
- return useTree;
-}
-
-void Clustal::promptForNewGuideTreeName(int type, string* treeName, const string& path)
-{
- string* ptrToMsg;
- if(type == Profile1)
- {
- ptrToMsg = &newProfile1TreePrompt;
- }
- else if(type == Profile2)
- {
- ptrToMsg = &newProfile2TreePrompt;
- }
- else
- {
- ptrToMsg = &newProfile1TreePrompt;
- }
-
- if(!treeName)
- {
- treeName = new string("");
- }
-
- while(treeName->empty())
- {
- string message = *ptrToMsg + path + "dnd]";
- string answer;
- utilityObject->getStr(message, answer);
- if(answer.empty())
- {
- answer = path + "dnd";
- *treeName = answer;
- }
- else
- {
- *treeName = answer;
- }
- }
-}
-
-}