clustercommand.cpp

   1 /*
   2  *  clustercommand.cpp
   3  *  Dotur
   4  *
   5  *  Created by Sarah Westcott on 1/2/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "clustercommand.h"
  11
  12 //**********************************************************************************************************************
  13 //This function checks to make sure the cluster command has no errors and then clusters based on the method chosen.
  14 ClusterCommand::ClusterCommand(string option){
  15         try{
  16                 globaldata = GlobalData::getInstance();
  17                 abort = false;
  18
  19                 //allow user to run help
  20                 if(option == "help") { help(); abort = true; }
  21
  22                 else {
  23                         //valid paramters for this command
  24                         string Array[] =  {"cutoff","precision","method"};
  25                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  26
  27                         OptionParser parser(option);
  28                         map<string,string> parameters = parser.getParameters();
  29
  30                         ValidParameters validParameter;
  31
  32                         //check to make sure all parameters are valid for command
  33                         for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
  34                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  35                         }
  36
  37                         //error checking to make sure they read a distance file
  38                         if ((globaldata->gSparseMatrix == NULL) || (globaldata->gListVector == NULL)) {
  39                                 mothurOut("Before you use the cluster command, you first need to read in a distance matrix."); mothurOutEndLine(); abort = true;
  40                         }
  41
  42                         //check for optional parameter and set defaults
  43                         // ...at some point should added some additional type checking...
  44                         //get user cutoff and precision or use defaults
  45                         string temp;
  46                         temp = validParameter.validFile(parameters, "precision", false);                if (temp == "not found") { temp = "100"; }
  47                         //saves precision legnth for formatting below
  48                         length = temp.length();
  49                         convert(temp, precision);
  50
  51                         temp = validParameter.validFile(parameters, "cutoff", false);                   if (temp == "not found") { temp = "10"; }
  52                         convert(temp, cutoff);
  53                         cutoff += (5 / (precision * 10.0));
  54
  55                         method = validParameter.validFile(parameters, "method", false);                 if (method == "not found") { method = "furthest"; }
  56
  57
  58                         if ((method == "furthest") || (method == "nearest") || (method == "average")) { }
  59                         else { mothurOut("Not a valid clustering method.  Valid clustering algorithms are furthest, nearest or average."); mothurOutEndLine(); abort = true; }
  60
  61
  62                         if (abort == false) {
  63
  64                                 //get matrix, list and rabund for execute
  65                                 if(globaldata->gSparseMatrix != NULL)   {       matrix = globaldata->gSparseMatrix;             }
  66
  67                                 if(globaldata->gListVector != NULL){
  68                                         list = globaldata->gListVector;
  69                                         rabund = new RAbundVector(list->getRAbundVector());
  70                                 }
  71
  72                                 //create cluster
  73                                 if(method == "furthest")        {       cluster = new CompleteLinkage(rabund, list, matrix);    tag = "fn";     }
  74                                 else if(method == "nearest"){   cluster = new SingleLinkage(rabund, list, matrix);              tag = "nn";     }
  75                                 else if(method == "average"){   cluster = new AverageLinkage(rabund, list, matrix);             tag = "an";     }
  76                                 else                                            {       mothurOut("error - not recognized method"); mothurOutEndLine(); abort = true;   }
  77
  78                                 fileroot = getRootName(globaldata->inputFileName);
  79
  80                                 openOutputFile(fileroot+ tag + ".sabund",       sabundFile);
  81                                 openOutputFile(fileroot+ tag + ".rabund",       rabundFile);
  82                                 openOutputFile(fileroot+ tag + ".list",         listFile);
  83
  84
  85                         }
  86
  87                 }
  88
  89         }
  90         catch(exception& e) {
  91                 errorOut(e, "ClusterCommand", "ClusterCommand");
  92                 exit(1);
  93         }
  94 }
  95
  96 //**********************************************************************************************************************
  97
  98 void ClusterCommand::help(){
  99         try {
 100                  mothurOut("The cluster command can only be executed after a successful read.dist command.\n");
 101                  mothurOut("The cluster command parameter options are method, cuttoff and precision. No parameters are required.\n");
 102                  mothurOut("The cluster command should be in the following format: \n");
 103                  mothurOut("cluster(method=yourMethod, cutoff=yourCutoff, precision=yourPrecision) \n");
 104                  mothurOut("The acceptable cluster methods are furthest, nearest and average.  If no method is provided then furthest is assumed.\n\n");
 105         }
 106         catch(exception& e) {
 107                 errorOut(e, "ClusterCommand", "help");
 108                 exit(1);
 109         }
 110 }
 111
 112 //**********************************************************************************************************************
 113
 114 ClusterCommand::~ClusterCommand(){
 115         if (abort == false) {
 116                 delete cluster;
 117                 delete rabund;
 118         }
 119 }
 120
 121 //**********************************************************************************************************************
 122
 123 int ClusterCommand::execute(){
 124         try {
 125
 126                 if (abort == true) {    return 0;       }
 127
 128                 float previousDist = 0.00000;
 129                 float rndPreviousDist = 0.00000;
 130                 oldRAbund = *rabund;
 131                 oldList = *list;
 132
 133                 float x;
 134                 x=0.1;
 135                 toString(x, 2);
 136
 137                 while(matrix->getSmallDist() < cutoff && matrix->getNNodes() > 0){
 138                         cluster->update();
 139                         float dist = matrix->getSmallDist();
 140                         float rndDist = roundDist(dist, precision);
 141
 142                         if(previousDist <= 0.0000 && dist != previousDist){
 143                                 printData("unique");
 144                         }
 145                         else if(rndDist != rndPreviousDist){
 146                                 printData(toString(rndPreviousDist,  length-1));
 147                         }
 148
 149                         previousDist = dist;
 150                         rndPreviousDist = rndDist;
 151                         oldRAbund = *rabund;
 152                         oldList = *list;
 153                 }
 154
 155                 if(previousDist <= 0.0000){
 156                         printData("unique");
 157                 }
 158                 else if(rndPreviousDist<cutoff){
 159                         printData(toString(rndPreviousDist, length-1));
 160                 }
 161
 162                 //delete globaldata's copy of the sparsematrix and listvector to free up memory
 163                 delete globaldata->gSparseMatrix;  globaldata->gSparseMatrix = NULL;
 164                 delete globaldata->gListVector;  globaldata->gListVector = NULL;
 165
 166                 //saves .list file so you can do the collect, rarefaction and summary commands without doing a read.list
 167                 if (globaldata->getFormat() == "phylip") { globaldata->setPhylipFile(""); }
 168                 else if (globaldata->getFormat() == "column") { globaldata->setColumnFile(""); }
 169
 170                 globaldata->setListFile(fileroot+ tag + ".list");
 171                 globaldata->setNameFile("");
 172                 globaldata->setFormat("list");
 173
 174                 sabundFile.close();
 175                 rabundFile.close();
 176                 listFile.close();
 177
 178                 return 0;
 179         }
 180         catch(exception& e) {
 181                 errorOut(e, "ClusterCommand", "execute");
 182                 exit(1);
 183         }
 184 }
 185
 186 //**********************************************************************************************************************
 187
 188 void ClusterCommand::printData(string label){
 189         try {
 190                 oldRAbund.setLabel(label);
 191                 oldRAbund.getSAbundVector().print(cout);
 192                 oldRAbund.print(rabundFile);
 193                 oldRAbund.getSAbundVector().print(sabundFile);
 194
 195                 oldList.setLabel(label);
 196                 oldList.print(listFile);
 197         }
 198         catch(exception& e) {
 199                 errorOut(e, "ClusterCommand", "printData");
 200                 exit(1);
 201         }
 202 }
 203 //**********************************************************************************************************************