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(){
  20         try {
  21                 globaldata = GlobalData::getInstance();
  22                 validCalculator = new ValidCalculators();
  23                 countends = globaldata->getCountEnds();
  24                 convert(globaldata->getProcessors(), processors);
  25                 convert(globaldata->getCutOff(), cutoff);
  26
  27                 //open file
  28                 string filename = globaldata->getFastaFile();
  29                 openInputFile(filename, in);
  30
  31
  32                 int i;
  33                 if (isTrue(countends) == true) {
  34                         for (i=0; i<globaldata->Estimators.size(); i++) {
  35                                 if (validCalculator->isValidCalculator("distance", globaldata->Estimators[i]) == true) {
  36                                         if (globaldata->Estimators[i] == "nogaps") {
  37                                                 distCalculator = new ignoreGaps();
  38                                         }else if (globaldata->Estimators[i] == "eachgap") {
  39                                                 distCalculator = new eachGapDist();
  40                                         }else if (globaldata->Estimators[i] == "onegap") {
  41                                         distCalculator = new oneGapDist();                                      }
  42                                 }
  43                         }
  44                 }else {
  45                         for (i=0; i<globaldata->Estimators.size(); i++) {
  46                                 if (validCalculator->isValidCalculator("distance", globaldata->Estimators[i]) == true) {
  47                                         if (globaldata->Estimators[i] == "nogaps") {
  48                                                 distCalculator = new ignoreGaps();
  49                                         }else if (globaldata->Estimators[i] == "eachgap") {
  50                                                 distCalculator = new eachGapIgnoreTermGapDist();
  51                                         }else if (globaldata->Estimators[i] == "onegap") {
  52                                                 distCalculator = new oneGapIgnoreTermGapDist();
  53                                         }
  54                                 }
  55                         }
  56                 }
  57
  58                 //reset calc for next command
  59                 globaldata->setCalc("");
  60         }
  61         catch(exception& e) {
  62                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function DistanceCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  63                 exit(1);
  64         }
  65         catch(...) {
  66                 cout << "An unknown error has occurred in the DistanceCommand class function DistanceCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  67                 exit(1);
  68         }
  69 }
  70 //**********************************************************************************************************************
  71
  72 int DistanceCommand::execute(){
  73         try {
  74
  75                 //reads fasta file and fills sequenceDB
  76                 if(globaldata->getFastaFile() != "") {  seqDB = new SequenceDB(in);  }
  77                 else { cout << "Error no fasta file." << endl; return 0; }
  78
  79                 int numSeqs = seqDB->getNumSeqs();
  80                 cutoff += 0.005;
  81
  82                 string distFile = getRootName(globaldata->getFastaFile()) + "dist";
  83
  84                 remove(distFile.c_str());
  85
  86                 //#     if defined (_WIN32)
  87                 //figure out how to implement the fork and wait commands in windows
  88                 //      driver(distCalculator, seqDB, 0, numSeqs, distFile, cutoff);
  89                 //#     endif
  90
  91 #if defined (__APPLE__) || (__MACH__)
  92                 if(processors == 1){
  93                         driver(distCalculator, seqDB, 0, numSeqs, distFile, cutoff);
  94                 }
  95                 else if(processors == 2){
  96
  97                         int pid = fork();
  98                         if(pid > 0){
  99                                 driver(distCalculator, seqDB, 0, (numSeqs/sqrt(2)), distFile + "tempa", cutoff);
 100                                 appendFiles((distFile+"tempa"), distFile);
 101                                 remove((distFile + "tempa").c_str());
 102                         }
 103                         else{
 104                                 driver(distCalculator, seqDB, (numSeqs/sqrt(2)), numSeqs, distFile + "tempb", cutoff);
 105                                 appendFiles((distFile+"tempb"), distFile);
 106                                 remove((distFile + "tempb").c_str());
 107                         }
 108                         wait(NULL);
 109
 110                 }
 111                 else if(processors == 3){
 112                         int pid1 = fork();
 113                         if(pid1 > 0){
 114                                 int pid2 = fork();
 115                                 if(pid2 > 0){
 116                                         driver(distCalculator, seqDB, 0, sqrt(3) * numSeqs / 3, distFile + "tempa", cutoff);
 117                                         appendFiles(distFile+"tempa", distFile);
 118                                         remove((distFile + "tempa").c_str());
 119                                 }
 120                                 else{
 121                                         driver(distCalculator, seqDB, sqrt(3) * numSeqs / 3, sqrt(6) * numSeqs / 3, distFile + "tempb", cutoff);
 122                                         appendFiles(distFile+"tempb", distFile);
 123                                         remove((distFile + "tempb").c_str());
 124                                 }
 125                                 wait(NULL);
 126                         }
 127                         else{
 128                                 driver(distCalculator, seqDB, sqrt(6) * numSeqs / 3, numSeqs, distFile + "tempc", cutoff);
 129                                 appendFiles(distFile+"tempc", distFile);
 130                                 remove((distFile + "tempc").c_str());
 131                         }
 132                         wait(NULL);
 133                 }
 134                 else if(processors == 4){
 135                         int pid1 = fork();
 136                         if(pid1 > 0){
 137                                 int pid2 = fork();
 138                                 if(pid2 > 0){
 139                                         driver(distCalculator, seqDB, 0, numSeqs / 2, distFile + "tempa", cutoff);
 140                                         appendFiles(distFile+"tempa", distFile);
 141                                         remove((distFile + "tempa").c_str());
 142                                 }
 143                                 else{
 144                                         driver(distCalculator, seqDB, numSeqs / 2, (numSeqs/sqrt(2)), distFile + "tempb", cutoff);
 145                                         appendFiles(distFile+"tempb", distFile);
 146                                         remove((distFile + "tempb").c_str());
 147                                 }
 148                                 wait(NULL);
 149                         }
 150                         else{
 151                                 int pid3 = fork();
 152                                 if(pid3 > 0){
 153                                         driver(distCalculator, seqDB, (numSeqs/sqrt(2)), (sqrt(3) * numSeqs / 2), distFile + "tempc", cutoff);
 154                                         appendFiles(distFile+"tempc", distFile);
 155                                         remove((distFile + "tempc").c_str());
 156                                 }
 157                                 else{
 158                                         driver(distCalculator, seqDB, (sqrt(3) * numSeqs / 2), numSeqs, distFile + "tempd", cutoff);
 159                                         appendFiles(distFile+"tempd", distFile);
 160                                         remove((distFile + "tempd").c_str());
 161                                 }
 162                                 wait(NULL);
 163                         }
 164                         wait(NULL);
 165                 }
 166                 wait(NULL);
 167 #elif (linux) || (__linux)
 168                 if(processors == 1){
 169                         driver(distCalculator, seqDB, 0, numSeqs, distFile, cutoff);
 170                 }
 171                 else if(processors == 2){
 172
 173                         int pid = fork();
 174                         if(pid > 0){
 175                                 driver(distCalculator, seqDB, 0, (numSeqs/sqrt(2)), distFile + "tempa", cutoff);
 176                                 appendFiles((distFile+"tempa"), distFile);
 177                                 remove((distFile + "tempa").c_str());
 178                         }
 179                         else{
 180                                 driver(distCalculator, seqDB, (numSeqs/sqrt(2)), numSeqs, distFile + "tempb", cutoff);
 181                                 appendFiles((distFile+"tempb"), distFile);
 182                                 remove((distFile + "tempb").c_str());
 183                         }
 184                         wait();
 185
 186                 }
 187                 else if(processors == 3){
 188                         int pid1 = fork();
 189                         if(pid1 > 0){
 190                                 int pid2 = fork();
 191                                 if(pid2 > 0){
 192                                         driver(distCalculator, seqDB, 0, sqrt(3) * numSeqs / 3, distFile + "tempa", cutoff);
 193                                         appendFiles(distFile+"tempa", distFile);
 194                                         remove((distFile + "tempa").c_str());
 195                                 }
 196                                 else{
 197                                         driver(distCalculator, seqDB, sqrt(3) * numSeqs / 3, sqrt(6) * numSeqs / 3, distFile + "tempb", cutoff);
 198                                         appendFiles(distFile+"tempb", distFile);
 199                                         remove((distFile + "tempb").c_str());
 200                                 }
 201                                 wait();
 202                         }
 203                         else{
 204                                 driver(distCalculator, seqDB, sqrt(6) * numSeqs / 3, numSeqs, distFile + "tempc", cutoff);
 205                                 appendFiles(distFile+"tempc", distFile);
 206                                 remove((distFile + "tempc").c_str());
 207                         }
 208                         wait();
 209                 }
 210                 else if(processors == 4){
 211                         int pid1 = fork();
 212                         if(pid1 > 0){
 213                                 int pid2 = fork();
 214                                 if(pid2 > 0){
 215                                         driver(distCalculator, seqDB, 0, numSeqs / 2, distFile + "tempa", cutoff);
 216                                         appendFiles(distFile+"tempa", distFile);
 217                                         remove((distFile + "tempa").c_str());
 218                                 }
 219                                 else{
 220                                         driver(distCalculator, seqDB, numSeqs / 2, (numSeqs/sqrt(2)), distFile + "tempb", cutoff);
 221                                         appendFiles(distFile+"tempb", distFile);
 222                                         remove((distFile + "tempb").c_str());
 223                                 }
 224                                 wait();
 225                         }
 226                         else{
 227                                 int pid3 = fork();
 228                                 if(pid3 > 0){
 229                                         driver(distCalculator, seqDB, (numSeqs/sqrt(2)), (sqrt(3) * numSeqs / 2), distFile + "tempc", cutoff);
 230                                         appendFiles(distFile+"tempc", distFile);
 231                                         remove((distFile + "tempc").c_str());
 232                                 }
 233                                 else{
 234                                         driver(distCalculator, seqDB, (sqrt(3) * numSeqs / 2), numSeqs, distFile + "tempd", cutoff);
 235                                         appendFiles(distFile+"tempd", distFile);
 236                                         remove((distFile + "tempd").c_str());
 237                                 }
 238                                 wait();
 239                         }
 240                         wait();
 241                 }
 242                 wait();
 243
 244 #else
 245                 driver(distCalculator, seqDB, 0, numSeqs, distFile, cutoff);
 246 #endif
 247
 248                 delete distCalculator;
 249
 250                 return 0;
 251
 252         }
 253         catch(exception& e) {
 254                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 255                 exit(1);
 256         }
 257         catch(...) {
 258                 cout << "An unknown error has occurred in the DistanceCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 259                 exit(1);
 260         }
 261 }
 262
 263 /**************************************************************************************************/
 264 /////// need to fix to work with calcs and sequencedb
 265 int DistanceCommand::driver(Dist* distCalculator, SequenceDB* align, int startLine, int endLine, string dFileName, float cutoff){
 266         try {
 267                 int startTime = time(NULL);
 268
 269                 ofstream distFile(dFileName.c_str(), ios::trunc);
 270                 distFile.setf(ios::fixed, ios::showpoint);
 271                 distFile << setprecision(4);
 272
 273                 for(int i=startLine;i<endLine;i++){
 274
 275                         for(int j=0;j<i;j++){
 276                                 distCalculator->calcDist(*(align->get(i)), *(align->get(j)));
 277                                 double dist = distCalculator->getDist();
 278
 279                                 if(dist <= cutoff){
 280                                         distFile << align->get(i)->getName() << ' ' << align->get(j)->getName() << ' ' << dist << endl;
 281                                 }
 282
 283                         }
 284                         if(i % 100 == 0){
 285                                 cout << i << '\t' << time(NULL) - startTime << endl;
 286                         }
 287
 288                 }
 289                 cout << endLine-1 << '\t' << time(NULL) - startTime << endl;
 290
 291                 return 1;
 292         }
 293         catch(exception& e) {
 294                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function driver. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 295                 exit(1);
 296         }
 297         catch(...) {
 298                 cout << "An unknown error has occurred in the DistanceCommand class function driver. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 299                 exit(1);
 300         }
 301
 302 }
 303
 304 /**************************************************************************************************/
 305 void DistanceCommand::appendFiles(string temp, string filename) {
 306         try{
 307                 ofstream output;
 308                 ifstream input;
 309
 310                 //open output file in append mode
 311                 openOutputFileAppend(filename, output);
 312
 313                 //open temp file for reading
 314                 openInputFile(temp, input);
 315
 316                 string line;
 317                 //read input file and write to output file
 318                 while(input.eof() != true) {
 319                         getline(input, line); //getline removes the newline char
 320                         if (line != "") {
 321                                 output << line << endl;   // Appending back newline char
 322                         }
 323                 }
 324
 325                 input.close();
 326                 output.close();
 327         }
 328         catch(exception& e) {
 329                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 330                 exit(1);
 331         }
 332         catch(...) {
 333                 cout << "An unknown error has occurred in the DistanceCommand class function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 334                 exit(1);
 335         }
 336 }
 337 /**************************************************************************************************/