binaries/src/clustalw/src/tree/UPGMA/RootedClusterTree.cpp

   1 /**
   2  * Author: Mark Larkin
   3  *
   4  * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
   5  */
   6 #ifdef HAVE_CONFIG_H
   7     #include "config.h"
   8 #endif
   9 #include "RootedClusterTree.h"
  10 #include "../../general/OutputFile.h"
  11 #include "UPGMAAlgorithm.h"
  12 #include "RootedTreeOutput.h"
  13 //#include "../RandomGenerator.h"
  14 namespace clustalw
  15 {
  16
  17 auto_ptr<AlignmentSteps> RootedClusterTree::treeFromDistMatrix(RootedGuideTree* phyloTree,DistMatrix* distMat, Alignment *alignPtr,
  18                                            int seq1, int nSeqs, string& phylipName)
  19 {
  20     OutputFile phylipPhyTreeFile;
  21     auto_ptr<AlignmentSteps> progSteps;
  22     try
  23     {
  24         // Test to see if the inputs are valid
  25         if(seq1 < 1 || nSeqs < 1)
  26         {
  27             cerr << "Invalid inputs into treeFromDistMatrix \n"
  28                  << "seq1 = " << seq1 << " nSeqs = " << nSeqs << "\n"
  29                  << "Need to end program!\n";
  30             exit(1);
  31             return progSteps;
  32         }
  33
  34         float dist;
  35         string path;
  36         verbose = false;
  37         firstSeq = seq1;
  38         lastSeq = firstSeq + nSeqs - 1;
  39
  40         SeqInfo info;
  41         info.firstSeq = firstSeq;
  42         info.lastSeq = lastSeq;
  43         info.numSeqs = nSeqs;
  44
  45         utilityObject->getPath(userParameters->getSeqName(), &path);
  46
  47         if(nSeqs >= 2)
  48         {
  49             string name = phylipName;
  50             if(!phylipPhyTreeFile.openFile(&name,
  51                              "\nEnter name for new GUIDE TREE           file  ", &path, "dnd",
  52                              "Guide tree"))
  53             {
  54                 return progSteps;
  55             }
  56             phylipName = name;
  57         }
  58         else
  59         {
  60             return progSteps;
  61         }
  62
  63         RootedTreeOutput outputTree(&info);
  64
  65         ofstream* ptrToFile = phylipPhyTreeFile.getPtrToFile();
  66
  67         if (nSeqs == 2)
  68         {
  69             dist = (*distMat)(firstSeq, firstSeq + 1) / 2.0;
  70             if(ptrToFile->is_open())
  71             {
  72                 (*ptrToFile) <<  "(" << alignPtr->getName(firstSeq) << ":"
  73                              << setprecision(5)
  74                              << dist << "," << alignPtr->getName(firstSeq + 1) << ":"
  75                              << setprecision(5) << dist <<");\n";
  76             }
  77             progSteps.reset(new AlignmentSteps);
  78             vector<int> groups;
  79             groups.resize(nSeqs + 1, 0);
  80             groups[1] = 1;
  81             groups[2] = 2;
  82         }
  83         else
  84         {
  85             UPGMAAlgorithm clusAlgorithm;
  86             progSteps = clusAlgorithm.generateTree(phyloTree, distMat, &info, false);
  87             outputTree.printPhylipTree(phyloTree, ptrToFile, alignPtr, distMat);
  88         }
  89         return progSteps;
  90     }
  91     catch(const exception &ex)
  92     {
  93         cerr << "ERROR: Error has occured in treeFromDistMatrix. "
  94              << "Need to terminate program.\n"
  95              << ex.what();
  96         exit(1);
  97     }
  98     catch(...)
  99     {
 100         cerr << "ERROR: Error has occured in treeFromDistMatrix. "
 101              << "Need to terminate program.\n";
 102         exit(1);
 103     }
 104 }
 105
 106 void RootedClusterTree::treeFromAlignment(TreeNames* treeNames, Alignment *alignPtr)
 107 {
 108     try
 109     {
 110         OutputFile phylipPhyTreeFile;
 111         OutputFile clustalPhyTreeFile;
 112         OutputFile distancesPhyTreeFile;
 113         OutputFile nexusPhyTreeFile;
 114         OutputFile pimFile;
 115
 116         RootedGuideTree phyloTree;
 117
 118         string path;
 119         int j;
 120         int overspill = 0;
 121         int totalDists;
 122         numSeqs = alignPtr->getNumSeqs(); // NOTE class variable
 123
 124         /**
 125          * Check if numSeqs is ok
 126          */
 127         if(!checkIfConditionsMet(numSeqs, 2))
 128         {
 129             return;
 130         }
 131
 132         firstSeq = 1;
 133         lastSeq = numSeqs;
 134
 135         // The SeqInfo struct is passed to reduce the number of parameters passed!
 136         SeqInfo info;
 137         info.firstSeq = firstSeq;
 138         info.lastSeq = lastSeq;
 139         info.numSeqs = numSeqs;
 140
 141         RootedTreeOutput outputTree(&info); // No bootstrap!
 142
 143         utilityObject->getPath(userParameters->getSeqName(), &path);
 144
 145         /**
 146          * Open the required output files.
 147          */
 148         if(!openFilesForTreeFromAlignment(&clustalPhyTreeFile, &phylipPhyTreeFile,
 149                         &distancesPhyTreeFile, &nexusPhyTreeFile, &pimFile, treeNames, &path))
 150         {
 151             return; // Problem opeing one of the files, cannot continue!
 152         }
 153
 154         int _lenFirstSeq = alignPtr->getSeqLength(firstSeq);
 155
 156         bootPositions.clear();
 157         bootPositions.resize(_lenFirstSeq + 2);
 158
 159         for (j = 1; j <= _lenFirstSeq; ++j)
 160         {
 161             bootPositions[j] = j;
 162         }
 163
 164         /**
 165          * Calculate quickDist and overspill
 166          */
 167         overspill = calcQuickDistMatForAll(clustalPhyTreeFile.getPtrToFile(),
 168                        phylipPhyTreeFile.getPtrToFile(), nexusPhyTreeFile.getPtrToFile(),
 169                        pimFile.getPtrToFile(), distancesPhyTreeFile.getPtrToFile(), alignPtr);
 170
 171         // check if any distances overflowed the distance corrections
 172         if (overspill > 0)
 173         {
 174             totalDists = (numSeqs *(numSeqs - 1)) / 2;
 175             overspillMessage(overspill, totalDists);
 176         }
 177
 178         if (userParameters->getOutputTreeClustal())
 179         {
 180             verbose = true;
 181         }
 182         // Turn on file output
 183
 184
 185         if (userParameters->getOutputTreeClustal() ||
 186             userParameters->getOutputTreePhylip()
 187             || userParameters->getOutputTreeNexus())
 188         {
 189             UPGMAAlgorithm clusAlgorithm;
 190             clusAlgorithm.setVerbose(true);
 191             clusAlgorithm.generateTree(&phyloTree, quickDistMat.get(), &info, false,
 192                                        clustalPhyTreeFile.getPtrToFile());
 193             clusAlgorithm.setVerbose(false);
 194         }
 195
 196         if (userParameters->getOutputTreePhylip())
 197         {
 198             outputTree.printPhylipTree(&phyloTree, phylipPhyTreeFile.getPtrToFile(), alignPtr,
 199                                         quickDistMat.get());
 200         }
 201
 202         if (userParameters->getOutputTreeNexus())
 203         {
 204             outputTree.printNexusTree(&phyloTree, nexusPhyTreeFile.getPtrToFile(), alignPtr,
 205                                        quickDistMat.get());
 206         }
 207
 208         /** Free up resources!!!!! */
 209
 210         treeGaps.clear();
 211         bootPositions.clear();
 212
 213     }
 214     catch(const exception& ex)
 215     {
 216         cerr << ex.what() << endl;
 217         utilityObject->error("Terminating program. Cannot continue\n");
 218         exit(1);
 219     }
 220 }
 221
 222 }