pairwiseseqscommand.cpp

   1 /*
   2  *  pairwiseseqscommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 10/20/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "pairwiseseqscommand.h"
  11
  12 //**********************************************************************************************************************
  13 vector<string> PairwiseSeqsCommand::setParameters(){
  14         try {
  15                 CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none","phylip-column",false,true,true); parameters.push_back(pfasta);
  16                 CommandParameter palign("align", "Multiple", "needleman-gotoh-blast-noalign", "needleman", "", "", "","",false,false); parameters.push_back(palign);
  17                 CommandParameter pmatch("match", "Number", "", "1.0", "", "", "","",false,false); parameters.push_back(pmatch);
  18                 CommandParameter pmismatch("mismatch", "Number", "", "-1.0", "", "", "","",false,false); parameters.push_back(pmismatch);
  19                 CommandParameter pgapopen("gapopen", "Number", "", "-2.0", "", "", "","",false,false); parameters.push_back(pgapopen);
  20                 CommandParameter pgapextend("gapextend", "Number", "", "-1.0", "", "", "","",false,false); parameters.push_back(pgapextend);
  21                 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
  22                 CommandParameter poutput("output", "Multiple", "column-lt-square-phylip", "column", "", "", "","phylip-column",false,false,true); parameters.push_back(poutput);
  23                 CommandParameter pcalc("calc", "Multiple", "nogaps-eachgap-onegap", "onegap", "", "", "","",false,false); parameters.push_back(pcalc);
  24                 CommandParameter pcountends("countends", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(pcountends);
  25                 CommandParameter pcompress("compress", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(pcompress);
  26                 CommandParameter pcutoff("cutoff", "Number", "", "1.0", "", "", "","",false,false,true); parameters.push_back(pcutoff);
  27                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
  28                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
  29
  30                 vector<string> myArray;
  31                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  32                 return myArray;
  33         }
  34         catch(exception& e) {
  35                 m->errorOut(e, "PairwiseSeqsCommand", "setParameters");
  36                 exit(1);
  37         }
  38 }
  39 //**********************************************************************************************************************
  40 string PairwiseSeqsCommand::getHelpString(){
  41         try {
  42                 string helpString = "";
  43                 helpString += "The pairwise.seqs command reads a fasta file and creates distance matrix.\n";
  44                 helpString += "The pairwise.seqs command parameters are fasta, align, match, mismatch, gapopen, gapextend, calc, output, cutoff and processors.\n";
  45                 helpString += "The fasta parameter is required. You may enter multiple fasta files by separating their names with dashes. ie. fasta=abrecovery.fasta-amzon.fasta \n";
  46                 helpString += "The align parameter allows you to specify the alignment method to use.  Your options are: gotoh, needleman, blast and noalign. The default is needleman.\n";
  47                 helpString += "The match parameter allows you to specify the bonus for having the same base. The default is 1.0.\n";
  48                 helpString += "The mistmatch parameter allows you to specify the penalty for having different bases.  The default is -1.0.\n";
  49                 helpString += "The gapopen parameter allows you to specify the penalty for opening a gap in an alignment. The default is -2.0.\n";
  50                 helpString += "The gapextend parameter allows you to specify the penalty for extending a gap in an alignment.  The default is -1.0.\n";
  51                 helpString += "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";
  52                 helpString += "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";
  53                 helpString += "The cutoff parameter allows you to specify maximum distance to keep. The default is 1.0.\n";
  54                 helpString += "The output parameter allows you to specify format of your distance matrix. Options are column, lt, and square. The default is column.\n";
  55                 helpString += "The compress parameter allows you to indicate that you want the resulting distance file compressed.  The default is false.\n";
  56                 helpString += "The pairwise.seqs command should be in the following format: \n";
  57                 helpString += "pairwise.seqs(fasta=yourfastaFile, align=yourAlignmentMethod) \n";
  58                 helpString += "Example pairwise.seqs(fasta=candidate.fasta, align=blast)\n";
  59                 helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFastaFile).\n";
  60                 return helpString;
  61         }
  62         catch(exception& e) {
  63                 m->errorOut(e, "PairwiseSeqsCommand", "getHelpString");
  64                 exit(1);
  65         }
  66 }
  67 //**********************************************************************************************************************
  68 string PairwiseSeqsCommand::getOutputPattern(string type) {
  69     try {
  70         string pattern = "";
  71
  72         if (type == "phylip") {  pattern = "[filename],[outputtag],dist"; }
  73         else if (type == "column") { pattern = "[filename],dist"; }
  74         else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
  75
  76         return pattern;
  77     }
  78     catch(exception& e) {
  79         m->errorOut(e, "PairwiseSeqsCommand", "getOutputPattern");
  80         exit(1);
  81     }
  82 }
  83 //**********************************************************************************************************************
  84 PairwiseSeqsCommand::PairwiseSeqsCommand(){
  85         try {
  86                 abort = true; calledHelp = true;
  87                 setParameters();
  88                 vector<string> tempOutNames;
  89                 outputTypes["phylip"] = tempOutNames;
  90                 outputTypes["column"] = tempOutNames;
  91         }
  92         catch(exception& e) {
  93                 m->errorOut(e, "PairwiseSeqsCommand", "PairwiseSeqsCommand");
  94                 exit(1);
  95         }
  96 }
  97 //**********************************************************************************************************************
  98 PairwiseSeqsCommand::PairwiseSeqsCommand(string option)  {
  99         try {
 100                 abort = false; calledHelp = false;
 101
 102                 //allow user to run help
 103                 if(option == "help") { help(); abort = true; calledHelp = true; }
 104                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
 105
 106                 else {
 107                         vector<string> myArray = setParameters();
 108
 109                         OptionParser parser(option);
 110                         map<string, string> parameters = parser.getParameters();
 111
 112                         ValidParameters validParameter("pairwise.seqs");
 113                         map<string, string>::iterator it;
 114
 115                         //check to make sure all parameters are valid for command
 116                         for (it = parameters.begin(); it != parameters.end(); it++) {
 117                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 118                         }
 119
 120                         //initialize outputTypes
 121                         vector<string> tempOutNames;
 122                         outputTypes["phylip"] = tempOutNames;
 123                         outputTypes["column"] = tempOutNames;
 124
 125                         //if the user changes the output directory command factory will send this info to us in the output parameter
 126                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = "";         }
 127
 128
 129                         //if the user changes the input directory command factory will send this info to us in the output parameter
 130                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 131
 132                         if (inputDir == "not found"){   inputDir = "";          }
 133
 134                         fastaFileName = validParameter.validFile(parameters, "fasta", false);
 135                         if (fastaFileName == "not found") {
 136                                 //if there is a current fasta file, use it
 137                                 string filename = m->getFastaFile();
 138                                 if (filename != "") { fastaFileNames.push_back(filename); m->mothurOut("Using " + filename + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
 139                                 else {  m->mothurOut("You have no current fastafile and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
 140                         }
 141                         else {
 142                                 m->splitAtDash(fastaFileName, fastaFileNames);
 143
 144                                 //go through files and make sure they are good, if not, then disregard them
 145                                 for (int i = 0; i < fastaFileNames.size(); i++) {
 146
 147                                         bool ignore = false;
 148                                         if (fastaFileNames[i] == "current") {
 149                                                 fastaFileNames[i] = m->getFastaFile();
 150                                                 if (fastaFileNames[i] != "") {  m->mothurOut("Using " + fastaFileNames[i] + " as input file for the fasta parameter where you had given current."); m->mothurOutEndLine(); }
 151                                                 else {
 152                                                         m->mothurOut("You have no current fastafile, ignoring current."); m->mothurOutEndLine(); ignore=true;
 153                                                         //erase from file list
 154                                                         fastaFileNames.erase(fastaFileNames.begin()+i);
 155                                                         i--;
 156                                                 }
 157                                         }
 158
 159                                         if (!ignore) {
 160
 161                                                 if (inputDir != "") {
 162                                                         string path = m->hasPath(fastaFileNames[i]);
 163                                                         //if the user has not given a path then, add inputdir. else leave path alone.
 164                                                         if (path == "") {       fastaFileNames[i] = inputDir + fastaFileNames[i];               }
 165                                                 }
 166
 167                                                 int ableToOpen;
 168                                                 ifstream in;
 169
 170                                                 ableToOpen = m->openInputFile(fastaFileNames[i], in, "noerror");
 171
 172                                                 //if you can't open it, try default location
 173                                                 if (ableToOpen == 1) {
 174                                                         if (m->getDefaultPath() != "") { //default path is set
 175                                                                 string tryPath = m->getDefaultPath() + m->getSimpleName(fastaFileNames[i]);
 176                                                                 m->mothurOut("Unable to open " + fastaFileNames[i] + ". Trying default " + tryPath); m->mothurOutEndLine();
 177                                                                 ifstream in2;
 178                                                                 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
 179                                                                 in2.close();
 180                                                                 fastaFileNames[i] = tryPath;
 181                                                         }
 182                                                 }
 183
 184                                                 //if you can't open it, try output location
 185                                                 if (ableToOpen == 1) {
 186                                                         if (m->getOutputDir() != "") { //default path is set
 187                                                                 string tryPath = m->getOutputDir() + m->getSimpleName(fastaFileNames[i]);
 188                                                                 m->mothurOut("Unable to open " + fastaFileNames[i] + ". Trying output directory " + tryPath); m->mothurOutEndLine();
 189                                                                 ifstream in2;
 190                                                                 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
 191                                                                 in2.close();
 192                                                                 fastaFileNames[i] = tryPath;
 193                                                         }
 194                                                 }
 195
 196                                                 in.close();
 197
 198                                                 if (ableToOpen == 1) {
 199                                                         m->mothurOut("Unable to open " + fastaFileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();
 200                                                         //erase from file list
 201                                                         fastaFileNames.erase(fastaFileNames.begin()+i);
 202                                                         i--;
 203                                                 }else {
 204                                                         m->setFastaFile(fastaFileNames[i]);
 205                                                 }
 206                                         }
 207                                 }
 208
 209                                 //make sure there is at least one valid file left
 210                                 if (fastaFileNames.size() == 0) { m->mothurOut("no valid files."); m->mothurOutEndLine(); abort = true; }
 211                         }
 212
 213                         //check for optional parameter and set defaults
 214                         // ...at some point should added some additional type checking...
 215                         string temp;
 216                         temp = validParameter.validFile(parameters, "match", false);            if (temp == "not found"){       temp = "1.0";                   }
 217                         m->mothurConvert(temp, match);
 218
 219                         temp = validParameter.validFile(parameters, "mismatch", false);         if (temp == "not found"){       temp = "-1.0";                  }
 220                         m->mothurConvert(temp, misMatch);
 221             if (misMatch > 0) { m->mothurOut("[ERROR]: mismatch must be negative.\n"); abort=true; }
 222
 223                         temp = validParameter.validFile(parameters, "gapopen", false);          if (temp == "not found"){       temp = "-2.0";                  }
 224                         m->mothurConvert(temp, gapOpen);
 225             if (gapOpen > 0) { m->mothurOut("[ERROR]: gapopen must be negative.\n"); abort=true; }
 226
 227                         temp = validParameter.validFile(parameters, "gapextend", false);        if (temp == "not found"){       temp = "-1.0";                  }
 228                         m->mothurConvert(temp, gapExtend);
 229             if (gapExtend > 0) { m->mothurOut("[ERROR]: gapextend must be negative.\n"); abort=true; }
 230
 231                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = m->getProcessors();      }
 232                         m->setProcessors(temp);
 233                         m->mothurConvert(temp, processors);
 234
 235                         temp = validParameter.validFile(parameters, "cutoff", false);           if(temp == "not found"){        temp = "1.0"; }
 236                         m->mothurConvert(temp, cutoff);
 237
 238                         temp = validParameter.validFile(parameters, "countends", false);        if(temp == "not found"){        temp = "T";     }
 239                         countends = m->isTrue(temp);
 240
 241                         temp = validParameter.validFile(parameters, "compress", false);         if(temp == "not found"){  temp = "F"; }
 242                         compress = m->isTrue(temp);
 243
 244                         align = validParameter.validFile(parameters, "align", false);           if (align == "not found"){      align = "needleman";    }
 245
 246                         output = validParameter.validFile(parameters, "output", false);         if(output == "not found"){      output = "column"; }
 247             if (output=="phylip") { output = "lt"; }
 248                         if ((output != "column") && (output != "lt") && (output != "square")) { m->mothurOut(output + " is not a valid output form. Options are column, lt and square. I will use column."); m->mothurOutEndLine(); output = "column"; }
 249
 250                         calc = validParameter.validFile(parameters, "calc", false);
 251                         if (calc == "not found") { calc = "onegap";  }
 252                         else {
 253                                  if (calc == "default")  {  calc = "onegap";  }
 254                         }
 255                         m->splitAtDash(calc, Estimators);
 256                 }
 257
 258         }
 259         catch(exception& e) {
 260                 m->errorOut(e, "PairwiseSeqsCommand", "PairwiseSeqsCommand");
 261                 exit(1);
 262         }
 263 }
 264 //**********************************************************************************************************************
 265
 266 int PairwiseSeqsCommand::execute(){
 267         try {
 268                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 269
 270                 longestBase = 2000; //will need to update this in driver if we find sequences with more bases.  hardcoded so we don't have the pre-read user fasta file.
 271
 272                 cutoff += 0.005;
 273
 274                 for (int s = 0; s < fastaFileNames.size(); s++) {
 275                         if (m->control_pressed) { outputTypes.clear(); return 0; }
 276
 277                         m->mothurOut("Processing sequences from " + fastaFileNames[s] + " ..." ); m->mothurOutEndLine();
 278
 279                         if (outputDir == "") {  outputDir += m->hasPath(fastaFileNames[s]); }
 280
 281                         ifstream inFASTA;
 282                         m->openInputFile(fastaFileNames[s], inFASTA);
 283                         alignDB = SequenceDB(inFASTA);
 284                         inFASTA.close();
 285
 286                         int numSeqs = alignDB.getNumSeqs();
 287                         int startTime = time(NULL);
 288                         string outputFile = "";
 289
 290             map<string, string> variables;
 291             variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s]));
 292                         if (output == "lt") { //does the user want lower triangle phylip formatted file
 293                                 variables["[outputtag]"] = "phylip";
 294                 outputFile = getOutputFileName("phylip", variables);
 295                                 m->mothurRemove(outputFile); outputTypes["phylip"].push_back(outputFile);
 296                         }else if (output == "column") { //user wants column format
 297                                 outputFile = getOutputFileName("column", variables);
 298                                 outputTypes["column"].push_back(outputFile);
 299                                 m->mothurRemove(outputFile);
 300                         }else { //assume square
 301                 variables["[outputtag]"] = "square";
 302                 outputFile = getOutputFileName("phylip", variables);
 303                                 m->mothurRemove(outputFile);
 304                                 outputTypes["phylip"].push_back(outputFile);
 305                         }
 306
 307                         #ifdef USE_MPI
 308
 309                         int pid, start, end;
 310                         int tag = 2001;
 311
 312                         MPI_Status status;
 313                         MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
 314                         MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 315
 316                         //each process gets where it should start and stop in the file
 317                         if (output != "square") {
 318                                 start = int (sqrt(float(pid)/float(processors)) * numSeqs);
 319                                 end = int (sqrt(float(pid+1)/float(processors)) * numSeqs);
 320                         }else{
 321                                 start = int ((float(pid)/float(processors)) * numSeqs);
 322                                 end = int ((float(pid+1)/float(processors)) * numSeqs);
 323                         }
 324
 325                         if (output == "column") {
 326                                 MPI_File outMPI;
 327                                 int amode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
 328
 329                         char filename[1024];
 330                                 strcpy(filename, outputFile.c_str());
 331
 332                                 MPI_File_open(MPI_COMM_WORLD, filename, amode, MPI_INFO_NULL, &outMPI);
 333
 334                                 if (pid == 0) { //you are the root process
 335
 336                                         //do your part
 337                                         string outputMyPart;
 338
 339                                         driverMPI(start, end, outMPI, cutoff);
 340
 341                                         if (m->control_pressed) { outputTypes.clear(); MPI_File_close(&outMPI);  m->mothurRemove(outputFile); return 0; }
 342
 343                                         //wait on chidren
 344                                         for(int i = 1; i < processors; i++) {
 345                                                 if (m->control_pressed) { outputTypes.clear();  MPI_File_close(&outMPI);   m->mothurRemove(outputFile);  return 0; }
 346
 347                                                 char buf[5];
 348                                                 MPI_Recv(buf, 5, MPI_CHAR, i, tag, MPI_COMM_WORLD, &status);
 349                                         }
 350                                 }else { //you are a child process
 351                                         //do your part
 352                                         driverMPI(start, end, outMPI, cutoff);
 353
 354                                         if (m->control_pressed) { outputTypes.clear();  MPI_File_close(&outMPI);  m->mothurRemove(outputFile);  return 0; }
 355
 356                                         char buf[5];
 357                                         strcpy(buf, "done");
 358                                         //tell parent you are done.
 359                                         MPI_Send(buf, 5, MPI_CHAR, 0, tag, MPI_COMM_WORLD);
 360                                 }
 361
 362                                 MPI_File_close(&outMPI);
 363
 364                         }else { //lower triangle format
 365                                 if (pid == 0) { //you are the root process
 366
 367                                         //do your part
 368                                         string outputMyPart;
 369                                         unsigned long long mySize;
 370
 371                                         if (output != "square"){ driverMPI(start, end, outputFile, mySize); }
 372                                         else { driverMPI(start, end, outputFile, mySize, output); }
 373
 374                                         if (m->control_pressed) {  outputTypes.clear();   m->mothurRemove(outputFile);   return 0; }
 375
 376                                         int amode=MPI_MODE_APPEND|MPI_MODE_WRONLY|MPI_MODE_CREATE; //
 377                                         MPI_File outMPI;
 378                                         MPI_File inMPI;
 379
 380                                         char filename[1024];
 381                                         strcpy(filename, outputFile.c_str());
 382
 383                                         MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
 384
 385                                         //wait on chidren
 386                                         for(int b = 1; b < processors; b++) {
 387                                                 unsigned long long fileSize;
 388
 389                                                 if (m->control_pressed) { outputTypes.clear();  MPI_File_close(&outMPI);  m->mothurRemove(outputFile);   return 0; }
 390
 391                                                 MPI_Recv(&fileSize, 1, MPI_LONG, b, tag, MPI_COMM_WORLD, &status);
 392
 393                                                 string outTemp = outputFile + toString(b) + ".temp";
 394
 395                                                 char* buf = new char[outTemp.length()];
 396                                                 memcpy(buf, outTemp.c_str(), outTemp.length());
 397
 398                                                 MPI_File_open(MPI_COMM_SELF, buf, MPI_MODE_DELETE_ON_CLOSE|MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);
 399                                                 delete buf;
 400
 401                                                 int count = 0;
 402                                                 while (count < fileSize) {
 403                                                         char buf2[1];
 404                                                         MPI_File_read(inMPI, buf2, 1, MPI_CHAR, &status);
 405                                                         MPI_File_write(outMPI, buf2, 1, MPI_CHAR, &status);
 406                                                         count += 1;
 407                                                 }
 408
 409                                                 MPI_File_close(&inMPI); //deleted on close
 410                                         }
 411
 412                                         MPI_File_close(&outMPI);
 413                                 }else { //you are a child process
 414                                         //do your part
 415                                         unsigned long long size;
 416                                         if (output != "square"){ driverMPI(start, end, (outputFile + toString(pid) + ".temp"), size); }
 417                                         else { driverMPI(start, end, (outputFile + toString(pid) + ".temp"), size, output); }
 418
 419                                         if (m->control_pressed) {  return 0; }
 420
 421                                         //tell parent you are done.
 422                                         MPI_Send(&size, 1, MPI_LONG, 0, tag, MPI_COMM_WORLD);
 423                                 }
 424                         }
 425                         MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
 426         #else
 427
 428                 //#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 429                         //if you don't need to fork anything
 430                         if(processors == 1){
 431                                 if (output != "square") {  driver(0, numSeqs, outputFile, cutoff); }
 432                                 else { driver(0, numSeqs, outputFile, "square");  }
 433                         }else{ //you have multiple processors
 434
 435                                 for (int i = 0; i < processors; i++) {
 436                                         distlinePair tempLine;
 437                                         lines.push_back(tempLine);
 438                                         if (output != "square") {
 439                                                 lines[i].start = int (sqrt(float(i)/float(processors)) * numSeqs);
 440                                                 lines[i].end = int (sqrt(float(i+1)/float(processors)) * numSeqs);
 441                                         }else{
 442                                                 lines[i].start = int ((float(i)/float(processors)) * numSeqs);
 443                                                 lines[i].end = int ((float(i+1)/float(processors)) * numSeqs);
 444                                         }
 445                                 }
 446
 447                                 createProcesses(outputFile);
 448                         }
 449                 //#else
 450                         //ifstream inFASTA;
 451                         //if (output != "square") {  driver(0, numSeqs, outputFile, cutoff); }
 452                         //else { driver(0, numSeqs, outputFile, "square");  }
 453                 //#endif
 454
 455         #endif
 456                         if (m->control_pressed) { outputTypes.clear();   m->mothurRemove(outputFile); return 0; }
 457
 458                         #ifdef USE_MPI
 459                                 MPI_Comm_rank(MPI_COMM_WORLD, &pid);
 460
 461                                 if (pid == 0) { //only one process should output to screen
 462                         #endif
 463
 464                         ifstream fileHandle;
 465                         fileHandle.open(outputFile.c_str());
 466                         if(fileHandle) {
 467                                 m->gobble(fileHandle);
 468                                 if (fileHandle.eof()) { m->mothurOut(outputFile + " is blank. This can result if there are no distances below your cutoff.");  m->mothurOutEndLine(); }
 469                         }
 470
 471                         if (compress) {
 472                                 m->mothurOut("Compressing..."); m->mothurOutEndLine();
 473                                 m->mothurOut("(Replacing " + outputFile + " with " + outputFile + ".gz)"); m->mothurOutEndLine();
 474                                 system(("gzip -v " + outputFile).c_str());
 475                                 outputNames.push_back(outputFile + ".gz");
 476                         }else { outputNames.push_back(outputFile); }
 477
 478                         #ifdef USE_MPI
 479                                 }
 480                         #endif
 481
 482                         m->mothurOut("It took " + toString(time(NULL) - startTime) + " to calculate the distances for " + toString(numSeqs) + " sequences."); m->mothurOutEndLine();
 483
 484                         if (m->control_pressed) { outputTypes.clear(); m->mothurRemove(outputFile); return 0; }
 485                 }
 486
 487                 //set phylip file as new current phylipfile
 488                 string current = "";
 489                 itTypes = outputTypes.find("phylip");
 490                 if (itTypes != outputTypes.end()) {
 491                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setPhylipFile(current); }
 492                 }
 493
 494                 //set column file as new current columnfile
 495                 itTypes = outputTypes.find("column");
 496                 if (itTypes != outputTypes.end()) {
 497                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setColumnFile(current); }
 498                 }
 499
 500                 m->mothurOutEndLine();
 501                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 502                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 503                 m->mothurOutEndLine();
 504
 505
 506                 return 0;
 507         }
 508         catch(exception& e) {
 509                 m->errorOut(e, "PairwiseSeqsCommand", "execute");
 510                 exit(1);
 511         }
 512 }
 513
 514 /**************************************************************************************************/
 515 void PairwiseSeqsCommand::createProcesses(string filename) {
 516         try {
 517         int process = 1;
 518                 processIDS.clear();
 519
 520 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 521
 522
 523                 //loop through and create all the processes you want
 524                 while (process != processors) {
 525                         int pid = fork();
 526
 527                         if (pid > 0) {
 528                                 processIDS.push_back(pid);
 529                                 process++;
 530                         }else if (pid == 0){
 531                                 if (output != "square") {  driver(lines[process].start, lines[process].end, filename + toString(getpid()) + ".temp", cutoff); }
 532                                 else { driver(lines[process].start, lines[process].end, filename + toString(getpid()) + ".temp", "square"); }
 533                                 exit(0);
 534                         }else {
 535                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
 536                                 for (int i=0;i<processIDS.size();i++) { int temp = processIDS[i]; kill (temp, SIGINT); }
 537                                 exit(0);
 538                         }
 539                 }
 540
 541                 //parent do my part
 542                 if (output != "square") {  driver(lines[0].start, lines[0].end, filename, cutoff); }
 543                 else { driver(lines[0].start, lines[0].end, filename, "square"); }
 544
 545
 546                 //force parent to wait until all the processes are done
 547                 for (int i=0;i<processIDS.size();i++) {
 548                         int temp = processIDS[i];
 549                         wait(&temp);
 550                 }
 551 #else
 552         //////////////////////////////////////////////////////////////////////////////////////////////////////
 553                 //Windows version shared memory, so be careful when passing variables through the distanceData struct.
 554                 //Above fork() will clone, so memory is separate, but that's not the case with windows,
 555                 //that's why the distance calculator was moved inside of the driver to make separate copies.
 556                 //////////////////////////////////////////////////////////////////////////////////////////////////////
 557
 558                 vector<pairwiseData*> pDataArray; //[processors-1];
 559                 DWORD   dwThreadIdArray[processors-1];
 560                 HANDLE  hThreadArray[processors-1];
 561
 562                 //Create processor-1 worker threads.
 563                 for( int i=0; i<processors-1; i++ ){
 564                         string extension = toString(i) + ".temp";
 565
 566                         // Allocate memory for thread data.
 567                         pairwiseData* tempDist = new pairwiseData((filename+extension), align, "square", Estimators[0], countends, output, alignDB, m, lines[i+1].start, lines[i+1].end, match, misMatch, gapOpen, gapExtend, longestBase, cutoff, i);
 568                         pDataArray.push_back(tempDist);
 569                         processIDS.push_back(i);
 570
 571                         if (output != "square") { hThreadArray[i] = CreateThread(NULL, 0, MyPairwiseThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);  }
 572             else { hThreadArray[i] = CreateThread(NULL, 0, MyPairwiseSquareThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);  }
 573                 }
 574
 575                 //do your part
 576                 if (output != "square") {  driver(lines[0].start, lines[0].end, filename, cutoff); }
 577                 else { driver(lines[0].start, lines[0].end, filename, "square"); }
 578
 579                 //Wait until all threads have terminated.
 580                 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
 581
 582                 //Close all thread handles and free memory allocations.
 583                 for(int i=0; i < pDataArray.size(); i++){
 584             if (pDataArray[i]->count != (pDataArray[i]->end-pDataArray[i]->start)) {
 585                 m->mothurOut("[ERROR]: process " + toString(i) + " only processed " + toString(pDataArray[i]->count) + " of " + toString(pDataArray[i]->end-pDataArray[i]->start) + " sequences assigned to it, quitting. \n"); m->control_pressed = true;
 586             }
 587                         CloseHandle(hThreadArray[i]);
 588                         delete pDataArray[i];
 589                 }
 590
 591 #endif
 592
 593         //append and remove temp files
 594                 for (int i=0;i<processIDS.size();i++) {
 595                         m->appendFiles((filename + toString(processIDS[i]) + ".temp"), filename);
 596                         m->mothurRemove((filename + toString(processIDS[i]) + ".temp"));
 597                 }
 598
 599         }
 600         catch(exception& e) {
 601                 m->errorOut(e, "PairwiseSeqsCommand", "createProcesses");
 602                 exit(1);
 603         }
 604 }
 605
 606 /**************************************************************************************************/
 607 /////// need to fix to work with calcs and sequencedb
 608 int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, float cutoff){
 609         try {
 610
 611                 int startTime = time(NULL);
 612
 613         Alignment* alignment;
 614         if(align == "gotoh")                    {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase);                 }
 615                 else if(align == "needleman")   {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);                                }
 616                 else if(align == "blast")               {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
 617                 else if(align == "noalign")             {       alignment = new NoAlign();                                                                                                      }
 618                 else {
 619                         m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
 620                         m->mothurOutEndLine();
 621                         alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
 622                 }
 623
 624         ValidCalculators validCalculator;
 625         Dist* distCalculator;
 626         if (countends) {
 627             if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
 628                 if (Estimators[0] == "nogaps")                  {       distCalculator = new ignoreGaps();      }
 629                 else if (Estimators[0] == "eachgap")    {       distCalculator = new eachGapDist();     }
 630                 else if (Estimators[0] == "onegap")             {       distCalculator = new oneGapDist();      }
 631             }
 632         }else {
 633             if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
 634                 if (Estimators[0] == "nogaps")          {       distCalculator = new ignoreGaps();                                      }
 635                 else if (Estimators[0] == "eachgap"){   distCalculator = new eachGapIgnoreTermGapDist();        }
 636                 else if (Estimators[0] == "onegap")     {       distCalculator = new oneGapIgnoreTermGapDist();         }
 637             }
 638         }
 639
 640                 //column file
 641                 ofstream outFile(dFileName.c_str(), ios::trunc);
 642                 outFile.setf(ios::fixed, ios::showpoint);
 643                 outFile << setprecision(4);
 644
 645                 if((output == "lt") && startLine == 0){ outFile << alignDB.getNumSeqs() << endl;        }
 646
 647                 for(int i=startLine;i<endLine;i++){
 648                         if(output == "lt")      {
 649                                 string name = alignDB.get(i).getName();
 650                                 if (name.length() < 10) { //pad with spaces to make compatible
 651                                         while (name.length() < 10) {  name += " ";  }
 652                                 }
 653                                 outFile << name << '\t';
 654                         }
 655
 656                         for(int j=0;j<i;j++){
 657
 658                                 if (m->control_pressed) { outFile.close(); delete alignment; delete distCalculator; return 0;  }
 659
 660                                 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
 661                                         alignment->resize(alignDB.get(i).getUnaligned().length()+1);
 662                                 }
 663
 664                                 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
 665                                         alignment->resize(alignDB.get(j).getUnaligned().length()+1);
 666                                 }
 667
 668                                 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
 669                                 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
 670
 671                                 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
 672                                 seqI.setAligned(alignment->getSeqAAln());
 673                                 seqJ.setAligned(alignment->getSeqBAln());
 674
 675                                 distCalculator->calcDist(seqI, seqJ);
 676                                 double dist = distCalculator->getDist();
 677
 678                 if (m->debug) { m->mothurOut("[DEBUG]: " + seqI.getName() + '\t' +  alignment->getSeqAAln() + '\n' + seqJ.getName() + alignment->getSeqBAln() + '\n' + "distance = " + toString(dist) + "\n"); }
 679
 680                                 if(dist <= cutoff){
 681                                         if (output == "column") { outFile << alignDB.get(i).getName() << ' ' << alignDB.get(j).getName() << ' ' << dist << endl; }
 682                                 }
 683                                 if (output == "lt") {  outFile << dist << '\t'; }
 684                         }
 685
 686                         if (output == "lt") { outFile << endl; }
 687
 688                         if(i % 100 == 0){
 689                                 m->mothurOutJustToScreen(toString(i) + "\t" + toString(time(NULL) - startTime)+"\n");
 690                         }
 691
 692                 }
 693                 m->mothurOutJustToScreen(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)+"\n");
 694
 695                 outFile.close();
 696         delete alignment;
 697         delete distCalculator;
 698
 699                 return 1;
 700         }
 701         catch(exception& e) {
 702                 m->errorOut(e, "PairwiseSeqsCommand", "driver");
 703                 exit(1);
 704         }
 705 }
 706 /**************************************************************************************************/
 707 /////// need to fix to work with calcs and sequencedb
 708 int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, string square){
 709         try {
 710
 711                 int startTime = time(NULL);
 712
 713         Alignment* alignment;
 714         if(align == "gotoh")                    {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase);                 }
 715                 else if(align == "needleman")   {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);                                }
 716                 else if(align == "blast")               {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
 717                 else if(align == "noalign")             {       alignment = new NoAlign();                                                                                                      }
 718                 else {
 719                         m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
 720                         m->mothurOutEndLine();
 721                         alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
 722                 }
 723
 724         ValidCalculators validCalculator;
 725         Dist* distCalculator;
 726         if (countends) {
 727             if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
 728                 if (Estimators[0] == "nogaps")                  {       distCalculator = new ignoreGaps();      }
 729                 else if (Estimators[0] == "eachgap")    {       distCalculator = new eachGapDist();     }
 730                 else if (Estimators[0] == "onegap")             {       distCalculator = new oneGapDist();      }
 731             }
 732         }else {
 733             if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
 734                 if (Estimators[0] == "nogaps")          {       distCalculator = new ignoreGaps();                                      }
 735                 else if (Estimators[0] == "eachgap"){   distCalculator = new eachGapIgnoreTermGapDist();        }
 736                 else if (Estimators[0] == "onegap")     {       distCalculator = new oneGapIgnoreTermGapDist();         }
 737             }
 738         }
 739
 740                 //column file
 741                 ofstream outFile(dFileName.c_str(), ios::trunc);
 742                 outFile.setf(ios::fixed, ios::showpoint);
 743                 outFile << setprecision(4);
 744
 745                 if(startLine == 0){     outFile << alignDB.getNumSeqs() << endl;        }
 746
 747                 for(int i=startLine;i<endLine;i++){
 748
 749                         string name = alignDB.get(i).getName();
 750                         //pad with spaces to make compatible
 751                         if (name.length() < 10) { while (name.length() < 10) {  name += " ";  } }
 752
 753                         outFile << name << '\t';
 754
 755                         for(int j=0;j<alignDB.getNumSeqs();j++){
 756
 757                                 if (m->control_pressed) { outFile.close(); delete alignment; delete distCalculator; return 0;  }
 758
 759                                 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
 760                                         alignment->resize(alignDB.get(i).getUnaligned().length()+1);
 761                                 }
 762
 763                                 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
 764                                         alignment->resize(alignDB.get(j).getUnaligned().length()+1);
 765                                 }
 766
 767                                 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
 768                                 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
 769
 770                                 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
 771                                 seqI.setAligned(alignment->getSeqAAln());
 772                                 seqJ.setAligned(alignment->getSeqBAln());
 773
 774                                 distCalculator->calcDist(seqI, seqJ);
 775                                 double dist = distCalculator->getDist();
 776
 777                                 outFile << dist << '\t';
 778
 779                 if (m->debug) { m->mothurOut("[DEBUG]: " + seqI.getName() + '\t' +  alignment->getSeqAAln() + '\n' + seqJ.getName() + alignment->getSeqBAln() + '\n' + "distance = " + toString(dist) + "\n"); }
 780                         }
 781
 782                         outFile << endl;
 783
 784                         if(i % 100 == 0){
 785                                 m->mothurOutJustToScreen(toString(i) + "\t" + toString(time(NULL) - startTime)+"\n");
 786                         }
 787
 788                 }
 789                 m->mothurOutJustToScreen(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)+"\n");
 790
 791                 outFile.close();
 792         delete alignment;
 793         delete distCalculator;
 794
 795                 return 1;
 796         }
 797         catch(exception& e) {
 798                 m->errorOut(e, "PairwiseSeqsCommand", "driver");
 799                 exit(1);
 800         }
 801 }
 802 #ifdef USE_MPI
 803 /**************************************************************************************************/
 804 /////// need to fix to work with calcs and sequencedb
 805 int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, MPI_File& outMPI, float cutoff){
 806         try {
 807                 MPI_Status status;
 808                 int startTime = time(NULL);
 809
 810         Alignment* alignment;
 811         if(align == "gotoh")                    {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase);                 }
 812                 else if(align == "needleman")   {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);                                }
 813                 else if(align == "blast")               {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
 814                 else if(align == "noalign")             {       alignment = new NoAlign();                                                                                                      }
 815                 else {
 816                         m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
 817                         m->mothurOutEndLine();
 818                         alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
 819                 }
 820
 821         ValidCalculators validCalculator;
 822         Dist* distCalculator;
 823         if (countends) {
 824             if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
 825                 if (Estimators[0] == "nogaps")                  {       distCalculator = new ignoreGaps();      }
 826                 else if (Estimators[0] == "eachgap")    {       distCalculator = new eachGapDist();     }
 827                 else if (Estimators[0] == "onegap")             {       distCalculator = new oneGapDist();      }
 828             }
 829         }else {
 830             if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
 831                 if (Estimators[0] == "nogaps")          {       distCalculator = new ignoreGaps();                                      }
 832                 else if (Estimators[0] == "eachgap"){   distCalculator = new eachGapIgnoreTermGapDist();        }
 833                 else if (Estimators[0] == "onegap")     {       distCalculator = new oneGapIgnoreTermGapDist();         }
 834             }
 835         }
 836
 837                 string outputString = "";
 838
 839                 for(int i=startLine;i<endLine;i++){
 840
 841                         for(int j=0;j<i;j++){
 842
 843                                 if (m->control_pressed) { delete alignment; delete distCalculator; return 0;  }
 844
 845                                 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
 846                                         alignment->resize(alignDB.get(i).getUnaligned().length()+1);
 847                                 }
 848
 849                                 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
 850                                         alignment->resize(alignDB.get(j).getUnaligned().length()+1);
 851                                 }
 852
 853                                 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
 854                                 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
 855
 856                                 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
 857                                 seqI.setAligned(alignment->getSeqAAln());
 858                                 seqJ.setAligned(alignment->getSeqBAln());
 859
 860                                 distCalculator->calcDist(seqI, seqJ);
 861                                 double dist = distCalculator->getDist();
 862
 863                 if (m->debug) { cout << ("[DEBUG]: " + seqI.getName() + '\t' +  alignment->getSeqAAln() + '\n' + seqJ.getName() + alignment->getSeqBAln() + '\n' + "distance = " + toString(dist) + "\n"); }
 864
 865                                 if(dist <= cutoff){
 866                                          outputString += (alignDB.get(i).getName() + ' ' + alignDB.get(j).getName() + ' ' + toString(dist) + '\n');
 867                                 }
 868                         }
 869
 870                         if(i % 100 == 0){
 871                                 //m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
 872                                 cout << i << '\t' << (time(NULL) - startTime) << endl;
 873                         }
 874
 875
 876                         //send results to parent
 877                         int length = outputString.length();
 878
 879                         char* buf = new char[length];
 880                         memcpy(buf, outputString.c_str(), length);
 881
 882                         MPI_File_write_shared(outMPI, buf, length, MPI_CHAR, &status);
 883                         outputString = "";
 884                         delete buf;
 885
 886                 }
 887                 delete alignment;
 888         delete distCalculator;
 889                 return 1;
 890         }
 891         catch(exception& e) {
 892                 m->errorOut(e, "PairwiseSeqsCommand", "driverMPI");
 893                 exit(1);
 894         }
 895 }
 896 /**************************************************************************************************/
 897 /////// need to fix to work with calcs and sequencedb
 898 int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsigned long long& size){
 899         try {
 900                 MPI_Status status;
 901
 902                 MPI_File outMPI;
 903                 int amode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
 904
 905                 char filename[1024];
 906                 strcpy(filename, file.c_str());
 907
 908                 MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
 909
 910                 Alignment* alignment;
 911         if(align == "gotoh")                    {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase);                 }
 912                 else if(align == "needleman")   {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);                                }
 913                 else if(align == "blast")               {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
 914                 else if(align == "noalign")             {       alignment = new NoAlign();                                                                                                      }
 915                 else {
 916                         m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
 917                         m->mothurOutEndLine();
 918                         alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
 919                 }
 920
 921         ValidCalculators validCalculator;
 922         Dist* distCalculator;
 923         if (countends) {
 924             if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
 925                 if (Estimators[0] == "nogaps")                  {       distCalculator = new ignoreGaps();      }
 926                 else if (Estimators[0] == "eachgap")    {       distCalculator = new eachGapDist();     }
 927                 else if (Estimators[0] == "onegap")             {       distCalculator = new oneGapDist();      }
 928             }
 929         }else {
 930             if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
 931                 if (Estimators[0] == "nogaps")          {       distCalculator = new ignoreGaps();                                      }
 932                 else if (Estimators[0] == "eachgap"){   distCalculator = new eachGapIgnoreTermGapDist();        }
 933                 else if (Estimators[0] == "onegap")     {       distCalculator = new oneGapIgnoreTermGapDist();         }
 934             }
 935         }
 936
 937
 938                 string outputString = "";
 939                 size = 0;
 940
 941                 if(startLine == 0){     outputString += toString(alignDB.getNumSeqs()) + "\n";  }
 942
 943                 for(int i=startLine;i<endLine;i++){
 944
 945                         string name = alignDB.get(i).getName();
 946                         if (name.length() < 10) { //pad with spaces to make compatible
 947                                 while (name.length() < 10) {  name += " ";  }
 948                         }
 949                         outputString += name + "\t";
 950
 951                         for(int j=0;j<i;j++){
 952
 953                                 if (m->control_pressed) { delete alignment; delete distCalculator; return 0;  }
 954
 955                                 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
 956                                         alignment->resize(alignDB.get(i).getUnaligned().length()+1);
 957                                 }
 958
 959                                 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
 960                                         alignment->resize(alignDB.get(j).getUnaligned().length()+1);
 961                                 }
 962
 963                                 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
 964                                 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
 965
 966                                 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
 967                                 seqI.setAligned(alignment->getSeqAAln());
 968                                 seqJ.setAligned(alignment->getSeqBAln());
 969
 970                                 distCalculator->calcDist(seqI, seqJ);
 971                                 double dist = distCalculator->getDist();
 972
 973                 if (m->debug) { cout << ("[DEBUG]: " + seqI.getName() + '\t' +  alignment->getSeqAAln() + '\n' + seqJ.getName() + alignment->getSeqBAln() + '\n' + "distance = " + toString(dist) + "\n"); }
 974
 975                                 outputString += toString(dist) + "\t";
 976                         }
 977
 978                         outputString += "\n";
 979
 980                         //send results to parent
 981                         int length = outputString.length();
 982                         char* buf = new char[length];
 983                         memcpy(buf, outputString.c_str(), length);
 984
 985                         MPI_File_write(outMPI, buf, length, MPI_CHAR, &status);
 986                         size += outputString.length();
 987                         outputString = "";
 988                         delete buf;
 989                 }
 990
 991                 MPI_File_close(&outMPI);
 992         delete alignment;
 993         delete distCalculator;
 994
 995                 return 1;
 996         }
 997         catch(exception& e) {
 998                 m->errorOut(e, "PairwiseSeqsCommand", "driverMPI");
 999                 exit(1);
1000         }
1001 }
1002 /**************************************************************************************************/
1003 /////// need to fix to work with calcs and sequencedb
1004 int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsigned long long& size, string square){
1005         try {
1006                 MPI_Status status;
1007
1008                 MPI_File outMPI;
1009                 int amode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
1010
1011                 char filename[1024];
1012                 strcpy(filename, file.c_str());
1013
1014                 MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
1015
1016                 Alignment* alignment;
1017         if(align == "gotoh")                    {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase);                 }
1018                 else if(align == "needleman")   {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);                                }
1019                 else if(align == "blast")               {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
1020                 else if(align == "noalign")             {       alignment = new NoAlign();                                                                                                      }
1021                 else {
1022                         m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
1023                         m->mothurOutEndLine();
1024                         alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
1025                 }
1026
1027         ValidCalculators validCalculator;
1028         Dist* distCalculator;
1029         if (countends) {
1030             if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
1031                 if (Estimators[0] == "nogaps")                  {       distCalculator = new ignoreGaps();      }
1032                 else if (Estimators[0] == "eachgap")    {       distCalculator = new eachGapDist();     }
1033                 else if (Estimators[0] == "onegap")             {       distCalculator = new oneGapDist();      }
1034             }
1035         }else {
1036             if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
1037                 if (Estimators[0] == "nogaps")          {       distCalculator = new ignoreGaps();                                      }
1038                 else if (Estimators[0] == "eachgap"){   distCalculator = new eachGapIgnoreTermGapDist();        }
1039                 else if (Estimators[0] == "onegap")     {       distCalculator = new oneGapIgnoreTermGapDist();         }
1040             }
1041         }
1042
1043                 string outputString = "";
1044                 size = 0;
1045
1046                 if(startLine == 0){     outputString += toString(alignDB.getNumSeqs()) + "\n";  }
1047
1048                 for(int i=startLine;i<endLine;i++){
1049
1050                         string name = alignDB.get(i).getName();
1051                         if (name.length() < 10) { //pad with spaces to make compatible
1052                                 while (name.length() < 10) {  name += " ";  }
1053                         }
1054                         outputString += name + "\t";
1055
1056                         for(int j=0;j<alignDB.getNumSeqs();j++){
1057
1058                                 if (m->control_pressed) {  delete alignment; return 0;  }
1059
1060                                 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
1061                                         alignment->resize(alignDB.get(i).getUnaligned().length()+1);
1062                                 }
1063
1064                                 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
1065                                         alignment->resize(alignDB.get(j).getUnaligned().length()+1);
1066                                 }
1067
1068                                 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
1069                                 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
1070
1071                                 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
1072                                 seqI.setAligned(alignment->getSeqAAln());
1073                                 seqJ.setAligned(alignment->getSeqBAln());
1074
1075                                 distCalculator->calcDist(seqI, seqJ);
1076                                 double dist = distCalculator->getDist();
1077
1078                                 outputString += toString(dist) + "\t";
1079
1080                 if (m->debug) { cout << ("[DEBUG]: " + seqI.getName() + '\t' +  alignment->getSeqAAln() + '\n' + seqJ.getName() + alignment->getSeqBAln() + '\n' + "distance = " + toString(dist) + "\n"); }
1081                         }
1082
1083                         outputString += "\n";
1084
1085                         //send results to parent
1086                         int length = outputString.length();
1087                         char* buf = new char[length];
1088                         memcpy(buf, outputString.c_str(), length);
1089
1090                         MPI_File_write(outMPI, buf, length, MPI_CHAR, &status);
1091                         size += outputString.length();
1092                         outputString = "";
1093                         delete buf;
1094                 }
1095
1096                 MPI_File_close(&outMPI);
1097
1098         delete alignment;
1099                 return 1;
1100         }
1101         catch(exception& e) {
1102                 m->errorOut(e, "PairwiseSeqsCommand", "driverMPI");
1103                 exit(1);
1104         }
1105 }
1106 #endif
1107 /**************************************************************************************************/
1108