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