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                 ends = globaldata->getEnds();
  24                 convert(globaldata->getProcessors(), processors);
  25                 convert(globaldata->getCutOff(), cutoff);
  26
  27                 int i;
  28                 if (ends != "T") {
  29                         for (i=0; i<globaldata->Estimators.size(); i++) {
  30                                 if (validCalculator->isValidCalculator("distance", globaldata->Estimators[i]) == true) {
  31                                         if (globaldata->Estimators[i] == "nogaps") {
  32                                                 distCalculator = new ignoreGaps();
  33                                         }else if (globaldata->Estimators[i] == "eachgap") {
  34                                                 distCalculator = new eachGapDist();
  35                                         }else if (globaldata->Estimators[i] == "onegap") {
  36                                                 distCalculator = new oneGapDist();                                      }
  37                                 }
  38                         }
  39                 }else {
  40                         for (i=0; i<globaldata->Estimators.size(); i++) {
  41                                 if (validCalculator->isValidCalculator("distance", globaldata->Estimators[i]) == true) {
  42                                         if (globaldata->Estimators[i] == "nogaps") {
  43                                                 distCalculator = new ignoreGaps();
  44                                         }else if (globaldata->Estimators[i] == "eachgap") {
  45                                                 distCalculator = new eachGapIgnoreTermGapDist();
  46                                         }else if (globaldata->Estimators[i] == "onegap") {
  47                                                 distCalculator = new oneGapIgnoreTermGapDist();
  48                                         }
  49                                 }
  50                         }
  51                 }
  52
  53                 //reset calc for next command
  54                 globaldata->setCalc("");
  55         }
  56         catch(exception& e) {
  57                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function DistanceCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  58                 exit(1);
  59         }
  60         catch(...) {
  61                 cout << "An unknown error has occurred in the DistanceCommand class function DistanceCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  62                 exit(1);
  63         }
  64 }
  65 //**********************************************************************************************************************
  66
  67 int DistanceCommand::execute(){
  68         try {
  69
  70                 //read file
  71                 string filename = globaldata->inputFileName;
  72
  73                 if(globaldata->getFastaFile() != "") {
  74                         readSeqs =  new ReadFasta(filename); }
  75                 else if(globaldata->getNexusFile() != "") {
  76                         readSeqs = new ReadNexus(filename); }
  77                 else if(globaldata->getClustalFile() != "") {
  78                         readSeqs = new ReadClustal(filename); }
  79                 else if(globaldata->getPhylipFile() != "") {
  80                         readSeqs = new ReadPhylip(filename); }
  81
  82                 readSeqs->read();
  83                 seqDB = readSeqs->getDB();
  84
  85                 int numSeqs = seqDB->getNumSeqs();
  86
  87                 string distFile = getRootName(globaldata->getFastaFile()) + "dist";
  88
  89                 remove(distFile.c_str());
  90
  91         //      #       if defined (WIN_VERSION)
  92                 //      driver(distCalculator, seqDB, 0, numSeqs, distFile, cutoff);
  93         //      #       endif
  94
  95         //      #       if defined (LINUX_VERSION)
  96                         if(processors == 1){
  97                                 driver(distCalculator, seqDB, 0, numSeqs, distFile, cutoff);
  98                         }
  99                         else if(processors == 2){
 100
 101                                 int pid = fork();
 102                                 if(pid > 0){
 103                                         driver(distCalculator, seqDB, 0, (numSeqs/sqrt(2)), distFile + "tempa", cutoff);
 104                                         appendFiles((distFile+"tempa"), distFile);
 105                                         remove((distFile + "tempa").c_str());
 106                                 }
 107                                 else{
 108                                         driver(distCalculator, seqDB, (numSeqs/sqrt(2)), numSeqs, distFile + "tempb", cutoff);
 109                                         appendFiles((distFile+"tempb"), distFile);
 110                                         remove((distFile + "tempb").c_str());
 111                                 }
 112                                 wait(NULL);
 113
 114                         }
 115                         else if(processors == 3){
 116                                 int pid1 = fork();
 117                                 if(pid1 > 0){
 118                                         int pid2 = fork();
 119                                         if(pid2 > 0){
 120                                                 driver(distCalculator, seqDB, 0, sqrt(3) * numSeqs / 3, distFile + "tempa", cutoff);
 121                                                 appendFiles(distFile+"tempa", distFile);
 122                                                 remove((distFile + "tempa").c_str());
 123                                         }
 124                                         else{
 125                                                 driver(distCalculator, seqDB, sqrt(3) * numSeqs / 3, sqrt(6) * numSeqs / 3, distFile + "tempb", cutoff);
 126                                                 appendFiles(distFile+"tempb", distFile);
 127                                                 remove((distFile + "tempb").c_str());
 128                                         }
 129                                         wait(NULL);
 130                                 }
 131                                 else{
 132                                         driver(distCalculator, seqDB, sqrt(6) * numSeqs / 3, numSeqs, distFile + "tempc", cutoff);
 133                                         appendFiles(distFile+"tempc", distFile);
 134                                         remove((distFile + "tempc").c_str());
 135                                 }
 136                                 wait(NULL);
 137                         }
 138                         else if(processors == 4){
 139                                 int pid1 = fork();
 140                                 if(pid1 > 0){
 141                                         int pid2 = fork();
 142                                         if(pid2 > 0){
 143                                                 driver(distCalculator, seqDB, 0, numSeqs / 2, distFile + "tempa", cutoff);
 144                                                 appendFiles(distFile+"tempa", distFile);
 145                                                 remove((distFile + "tempa").c_str());
 146                                         }
 147                                         else{
 148                                                 driver(distCalculator, seqDB, numSeqs / 2, (numSeqs/sqrt(2)), distFile + "tempb", cutoff);
 149                                                 appendFiles(distFile+"tempb", distFile);
 150                                                 remove((distFile + "tempb").c_str());
 151                                         }
 152                                         wait(NULL);
 153                                 }
 154                                 else{
 155                                         int pid3 = fork();
 156                                         if(pid3 > 0){
 157                                                 driver(distCalculator, seqDB, (numSeqs/sqrt(2)), (sqrt(3) * numSeqs / 2), distFile + "tempc", cutoff);
 158                                                 appendFiles(distFile+"tempc", distFile);
 159                                                 remove((distFile + "tempc").c_str());
 160                                         }
 161                                         else{
 162                                                 driver(distCalculator, seqDB, (sqrt(3) * numSeqs / 2), numSeqs, distFile + "tempd", cutoff);
 163                                                 appendFiles(distFile+"tempd", distFile);
 164                                                 remove((distFile + "tempd").c_str());
 165                                         }
 166                                         wait(NULL);
 167                                 }
 168                                 wait(NULL);
 169                         }
 170                         wait(NULL);
 171                 //#     endif
 172
 173                 delete distCalculator;
 174
 175                 return 0;
 176
 177         }
 178         catch(exception& e) {
 179                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 180                 exit(1);
 181         }
 182         catch(...) {
 183                 cout << "An unknown error has occurred in the DistanceCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 184                 exit(1);
 185         }
 186 }
 187
 188 /**************************************************************************************************/
 189 /////// need to fix to work with calcs and sequencedb
 190 int DistanceCommand::driver(Dist* distCalculator, SequenceDB* align, int startLine, int endLine, string dFileName, float cutoff){
 191         try {
 192                 int startTime = time(NULL);
 193
 194                 ofstream distFile(dFileName.c_str(), ios::trunc);
 195                 distFile.setf(ios::fixed, ios::showpoint);
 196                 distFile << setprecision(4);
 197
 198                 for(int i=startLine;i<endLine;i++){
 199
 200                         for(int j=0;j<i;j++){
 201 //cout << "unaligned" << endl;
 202 //cout << align->get(i).getUnaligned() << "  " << align->get(j).getUnaligned() << endl;
 203 //cout << "aligned" << endl;
 204 //cout << align->get(i).getAligned() << "  " << align->get(j).getAligned() << endl;
 205                                 distCalculator->calcDist(align->get(i), align->get(j));
 206                                 double dist = distCalculator->getDist();
 207
 208                                 if(dist <= cutoff){
 209                                         distFile << align->get(i).getName() << ' ' << align->get(j).getName() << ' ' << dist << endl;
 210                                 }
 211
 212                         }
 213                         if(i % 100 == 0){
 214                                 cout << i << '\t' << time(NULL) - startTime << endl;
 215                         }
 216
 217                 }
 218                 cout << endLine-1 << '\t' << time(NULL) - startTime << endl;
 219
 220                 return 1;
 221         }
 222         catch(exception& e) {
 223                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function driver. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 224                 exit(1);
 225         }
 226         catch(...) {
 227                 cout << "An unknown error has occurred in the DistanceCommand class function driver. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 228                 exit(1);
 229         }
 230
 231 }
 232
 233 /**************************************************************************************************/
 234 void DistanceCommand::appendFiles(string temp, string filename) {
 235         try{
 236                 ofstream output;
 237                 ifstream input;
 238
 239                 //open output file in append mode
 240                 openOutputFileAppend(filename, output);
 241
 242                 //open temp file for reading
 243                 openInputFile(temp, input);
 244
 245                 string line;
 246                 //read input file and write to output file
 247                 while(input.eof() != true) {
 248                         getline(input, line); //getline removes the newline char
 249                         if (line != "") {
 250                                 output << line << endl;   // Appending back newline char
 251                         }
 252                 }
 253
 254                 input.close();
 255                 output.close();
 256         }
 257         catch(exception& e) {
 258                 cout << "Standard Error: " << e.what() << " has occurred in the DistanceCommand class Function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 259                 exit(1);
 260         }
 261         catch(...) {
 262                 cout << "An unknown error has occurred in the DistanceCommand class function appendFiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 263                 exit(1);
 264         }
 265 }
 266 /**************************************************************************************************/