distancecommand.cpp

   1 /*
   2  *  distancecommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 5/7/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "distancecommand.h"
  11 #include "ignoregaps.h"
  12 #include "eachgapdist.h"
  13 #include "eachgapignore.h"
  14 #include "onegapdist.h"
  15 #include "onegapignore.h"
  16
  17 //**********************************************************************************************************************
  18
  19 DistanceCommand::DistanceCommand(string option){
  20         try {
  21                 globaldata = GlobalData::getInstance();
  22                 abort = false;
  23                 Estimators.clear();
  24
  25                 //allow user to run help
  26                 if(option == "help") { help(); abort = true; }
  27
  28                 else {
  29                         //valid paramters for this command
  30                         string Array[] =  {"fasta", "phylip", "calc", "countends", "cutoff", "processors"};
  31                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  32
  33                         parser = new OptionParser();
  34                         parser->parse(option, parameters);  delete parser;
  35
  36                         ValidParameters* validParameter = new ValidParameters();
  37
  38                         //check to make sure all parameters are valid for command
  39                         for (it2 = parameters.begin(); it2 != parameters.end(); it2++) {
  40                                 if (validParameter->isValidParameter(it2->first, myArray, it2->second) != true) {  abort = true;  }
  41                         }
  42
  43                         //check for required parameters
  44                         fastafile = validParameter->validFile(parameters, "fasta", true);
  45                         if (fastafile == "not found") { cout << "fasta is a required parameter for the dist.seqs command." << endl; abort = true; }
  46                         else if (fastafile == "not open") { abort = true; }
  47                         else {
  48                                 globaldata->setFastaFile(fastafile);
  49                                 openInputFile(fastafile, in);
  50                         }
  51
  52                         //check for optional parameter and set defaults
  53                         // ...at some point should added some additional type checking...
  54                         calc = validParameter->validFile(parameters, "calc", false);
  55                         if (calc == "not found") { calc = "onegap";  }
  56                         else {
  57                                  if (calc == "default")  {  calc = "onegap";  }
  58                         }
  59                         splitAtDash(calc, Estimators);
  60
  61                         string temp;
  62                         temp = validParameter->validFile(parameters, "countends", false);                       if (temp == "not found") { temp = "T"; }
  63                         convert(temp, countends);
  64
  65                         temp = validParameter->validFile(parameters, "cutoff", false);                          if (temp == "not found") { temp = "1.0"; }
  66                         convert(temp, cutoff);
  67
  68                         temp = validParameter->validFile(parameters, "processors", false);                      if (temp == "not found") { temp = "1"; }
  69                         convert(temp, processors);
  70
  71                         phylip = validParameter->validFile(parameters, "phylip", false);                        if (phylip == "not found") { phylip = "F"; }
  72
  73                         delete validParameter;
  74
  75                         validCalculator = new ValidCalculators();
  76
  77                         int i;
  78                         if (isTrue(countends) == true) {
  79                                 for (i=0; i<Estimators.size(); i++) {
  80                                         if (validCalculator->isValidCalculator("distance", Estimators[i]) == true) {
  81                                                 if (Estimators[i] == "nogaps") {
  82                                                         distCalculator = new ignoreGaps();
  83                                                 }else if (Estimators[i] == "eachgap") {
  84                                                         distCalculator = new eachGapDist();
  85                                                 }else if (Estimators[i] == "onegap") {
  86                                                 distCalculator = new oneGapDist();                                      }
  87                                         }
  88                                 }
  89                         }else {
  90                                 for (i=0; i<Estimators.size(); i++) {
  91                                         if (validCalculator->isValidCalculator("distance", Estimators[i]) == true) {
  92                                                 if (Estimators[i] == "nogaps") {
  93                                                         distCalculator = new ignoreGaps();
  94                                                 }else if (Estimators[i] == "eachgap") {
  95                                                         distCalculator = new eachGapIgnoreTermGapDist();
  96                                                 }else if (Estimators[i] == "onegap") {
  97                                                         distCalculator = new oneGapIgnoreTermGapDist();
  98                                                 }
  99                                         }
 100                                 }
 101                         }
 102
 103                         delete validCalculator;
 104                 }
 105
 106         }
 107         catch(exception& e) {
 108                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function DistanceCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 109                 exit(1);
 110         }
 111         catch(...) {
 112                 cout << "An unknown error has occurred in the DistanceCommand class function DistanceCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 113                 exit(1);
 114         }
 115 }
 116 //**********************************************************************************************************************
 117
 118 void DistanceCommand::help(){
 119         try {
 120                 cout << "The dist.seqs command reads a file containing sequences and creates a distance file." << "\n";
 121                 cout << "The dist.seqs command parameters are fasta, calc, countends, cutoff and processors.  " << "\n";
 122                 cout << "The fasta parameter is required." << "\n";
 123                 cout << "The calc parameter allows you to specify the method of calculating the distances.  Your options are: nogaps, onegap or eachgap. The default is onegap." << "\n";
 124                 cout << "The countends parameter allows you to specify whether to include terminal gaps in distance.  Your options are: T or F. The default is T." << "\n";
 125                 cout << "The cutoff parameter allows you to specify maximum distance to keep. The default is 1.0." << "\n";
 126                 cout << "The processors parameter allows you to specify number of processors to use.  The default is 1." << "\n";
 127                 cout << "The dist.seqs command should be in the following format: " << "\n";
 128                 cout << "dist.seqs(fasta=yourFastaFile, calc=yourCalc, countends=yourEnds, cutoff= yourCutOff, processors=yourProcessors) " << "\n";
 129                 cout << "Example dist.seqs(fasta=amazon.fasta, calc=eachgap, countends=F, cutoff= 2.0, processors=3)." << "\n";
 130                 cout << "Note: No spaces between parameter labels (i.e. calc), '=' and parameters (i.e.yourCalc)." << "\n" << "\n";
 131         }
 132         catch(exception& e) {
 133                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function help. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 134                 exit(1);
 135         }
 136         catch(...) {
 137                 cout << "An unknown error has occurred in the DistanceCommand class function help. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 138                 exit(1);
 139         }
 140 }
 141
 142
 143 //**********************************************************************************************************************
 144
 145 int DistanceCommand::execute(){
 146         try {
 147
 148                 if (abort == true) { return 0; }
 149
 150                 //reads fasta file and fills sequenceDB
 151                 seqDB = new SequenceDB(in);
 152
 153                 int numSeqs = seqDB->getNumSeqs();
 154                 cutoff += 0.005;
 155
 156                 string outputFile;
 157
 158                 //doses the user want the phylip formatted file as well
 159                 if (isTrue(phylip) == true) {
 160                         outputFile = getRootName(fastafile) + "phylip.dist";
 161                         remove(outputFile.c_str());
 162
 163                         //output numSeqs to phylip formatted dist file
 164                         openOutputFile(outputFile, outFile);
 165                         outFile << numSeqs << endl;
 166                         outFile.close();
 167                 }else { //user wants column format
 168                         outputFile = getRootName(fastafile) + "dist";
 169                         remove(outputFile.c_str());
 170                 }
 171
 172                 //#     if defined (_WIN32)
 173                 //figure out how to implement the fork and wait commands in windows
 174                 //      driver(distCalculator, seqDB, 0, numSeqs, distFile, phylipFile, cutoff);
 175                 //#     endif
 176
 177
 178 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 179                 //if you don't need to fork anything
 180                 if(processors == 1){
 181                         driver(distCalculator, seqDB, 0, numSeqs, outputFile + ".temp", cutoff);
 182                         appendFiles((outputFile + ".temp"), outputFile);
 183                         remove((outputFile + ".temp").c_str());
 184                 }else{ //you have multiple processors
 185
 186                         for (int i = 0; i < processors; i++) {
 187                                 lines.push_back(new linePair());
 188                                 lines[i]->start = int (sqrt(float(i)/float(processors)) * numSeqs);
 189                                 lines[i]->end = int (sqrt(float(i+1)/float(processors)) * numSeqs);
 190                         }
 191
 192                         cout << lines[0]->start << '\t' << lines[0]->end << endl;
 193                         cout << lines[1]->start << '\t' << lines[1]->end << endl;
 194
 195                         createProcesses(outputFile);
 196
 197                         //append and remove temp files
 198                         for (it = processIDS.begin(); it != processIDS.end(); it++) {
 199                                 appendFiles((outputFile + toString(it->second) + ".temp"), outputFile);
 200                                 remove((outputFile + toString(it->second) + ".temp").c_str());
 201                         }
 202                 }
 203 #else
 204                 driver(distCalculator, seqDB, 0, numSeqs, outputFile + ".temp", cutoff);
 205                 appendFiles((outputFile + ".temp"), outputFile);
 206                 remove((outputFile + ".temp").c_str());
 207 #endif
 208
 209                 delete distCalculator;
 210
 211                 return 0;
 212
 213         }
 214         catch(exception& e) {
 215                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 216                 exit(1);
 217         }
 218         catch(...) {
 219                 cout << "An unknown error has occurred in the DistanceCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 220                 exit(1);
 221         }
 222 }
 223 /**************************************************************************************************/
 224 void DistanceCommand::createProcesses(string filename) {
 225         try {
 226                 int process = 0;
 227                 processIDS.clear();
 228
 229                 //loop through and create all the processes you want
 230                 while (process != processors) {
 231                         int pid = fork();
 232
 233                         if (pid > 0) {
 234                                 processIDS[lines[process]->end] = pid;  //create map from line number to pid so you can append files in correct order later
 235                                 process++;
 236                         }else if (pid == 0){
 237                                 driver(distCalculator, seqDB, lines[process]->start, lines[process]->end, filename + toString(getpid()) + ".temp", cutoff);
 238                                 exit(0);
 239                         }else { cout << "unable to spawn the necessary processes." << endl; exit(0); }
 240                 }
 241
 242                 //force parent to wait until all the processes are done
 243                 for (it = processIDS.begin(); it != processIDS.end(); it++) {
 244                         int temp = it->second;
 245                         wait(&temp);
 246                 }
 247
 248         }
 249         catch(exception& e) {
 250                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function createProcesses. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 251                 exit(1);
 252         }
 253         catch(...) {
 254                 cout << "An unknown error has occurred in the DistanceCommand class function createProcesses. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 255                 exit(1);
 256         }
 257 }
 258
 259 /**************************************************************************************************/
 260 /////// need to fix to work with calcs and sequencedb
 261 int DistanceCommand::driver(Dist* distCalculator, SequenceDB* align, int startLine, int endLine, string dFileName, float cutoff){
 262         try {
 263
 264                 int startTime = time(NULL);
 265
 266                 //column file
 267                 ofstream outFile(dFileName.c_str(), ios::trunc);
 268                 outFile.setf(ios::fixed, ios::showpoint);
 269                 outFile << setprecision(4);
 270
 271                 for(int i=startLine;i<endLine;i++){
 272
 273                         for(int j=0;j<i;j++){
 274                                 distCalculator->calcDist(*(align->get(i)), *(align->get(j)));
 275                                 double dist = distCalculator->getDist();
 276
 277                                 if(dist <= cutoff){
 278                                         if (isTrue(phylip) != true) { outFile << align->get(i)->getName() << ' ' << align->get(j)->getName() << ' ' << dist << endl; }
 279                                 }
 280                                 if (isTrue(phylip) == true) {  outFile << dist << '\t'; }
 281
 282                         }
 283
 284                         if (isTrue(phylip) == true) { outFile << endl; }
 285
 286                         if(i % 100 == 0){
 287                                 cout << i << '\t' << time(NULL) - startTime << endl;
 288                         }
 289
 290                 }
 291                 cout << endLine-1 << '\t' << time(NULL) - startTime << endl;
 292
 293                 //philFile.close();
 294                 //distFile.close();
 295
 296                 return 1;
 297         }
 298         catch(exception& e) {
 299                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function driver. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 300                 exit(1);
 301         }
 302         catch(...) {
 303                 cout << "An unknown error has occurred in the DistanceCommand class function driver. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 304                 exit(1);
 305         }
 306
 307 }
 308
 309 /**************************************************************************************************/
 310 void DistanceCommand::appendFiles(string temp, string filename) {
 311         try{
 312                 ofstream output;
 313                 ifstream input;
 314
 315                 //open output file in append mode
 316                 openOutputFileAppend(filename, output);
 317
 318                 //open temp file for reading
 319                 openInputFile(temp, input);
 320
 321                 string line;
 322                 //read input file and write to output file
 323                 while(input.eof() != true) {
 324                         getline(input, line); //getline removes the newline char
 325                         if (line != "") {
 326                                 output << line << endl;   // Appending back newline char
 327                         }
 328                 }
 329
 330                 input.close();
 331                 output.close();
 332         }
 333         catch(exception& e) {
 334                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 335                 exit(1);
 336         }
 337         catch(...) {
 338                 cout << "An unknown error has occurred in the DistanceCommand class function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 339                 exit(1);
 340         }
 341 }
 342 /**************************************************************************************************/