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