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