]> git.donarmstrong.com Git - mothur.git/blobdiff - pairwiseseqscommand.cpp
sffinfo bug with flow grams right index when clipQualRight=0
[mothur.git] / pairwiseseqscommand.cpp
index 1aa133f38c379fd40c495d2c2a844de292c17f51..c9e5ecfdfda727427793e98aaed0270cdcd6e046 100644 (file)
@@ -8,37 +8,24 @@
  */
 
 #include "pairwiseseqscommand.h"
-#include "sequence.hpp"
-
-#include "gotohoverlap.hpp"
-#include "needlemanoverlap.hpp"
-#include "blastalign.hpp"
-#include "noalign.hpp"
-
-#include "ignoregaps.h"
-#include "eachgapdist.h"
-#include "eachgapignore.h"
-#include "onegapdist.h"
-#include "onegapignore.h"
-
 
 //**********************************************************************************************************************
 vector<string> PairwiseSeqsCommand::setParameters(){   
        try {
-               CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pfasta);
-               CommandParameter palign("align", "Multiple", "needleman-gotoh-blast-noalign", "needleman", "", "", "",false,false); parameters.push_back(palign);
-               CommandParameter pmatch("match", "Number", "", "1.0", "", "", "",false,false); parameters.push_back(pmatch);
-               CommandParameter pmismatch("mismatch", "Number", "", "-1.0", "", "", "",false,false); parameters.push_back(pmismatch);
-               CommandParameter pgapopen("gapopen", "Number", "", "-2.0", "", "", "",false,false); parameters.push_back(pgapopen);
-               CommandParameter pgapextend("gapextend", "Number", "", "-1.0", "", "", "",false,false); parameters.push_back(pgapextend);
-               CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
-               CommandParameter poutput("output", "Multiple", "column-lt-square", "column", "", "", "",false,false); parameters.push_back(poutput);
-               CommandParameter pcalc("calc", "Multiple", "nogaps-eachgap-onegap", "onegap", "", "", "",false,false); parameters.push_back(pcalc);
-               CommandParameter pcountends("countends", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pcountends);
-               CommandParameter pcompress("compress", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pcompress);
-               CommandParameter pcutoff("cutoff", "Number", "", "1.0", "", "", "",false,false); parameters.push_back(pcutoff);
-               CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
-               CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
+               CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none","phylip-column",false,true,true); parameters.push_back(pfasta);
+               CommandParameter palign("align", "Multiple", "needleman-gotoh-blast-noalign", "needleman", "", "", "","",false,false); parameters.push_back(palign);
+               CommandParameter pmatch("match", "Number", "", "1.0", "", "", "","",false,false); parameters.push_back(pmatch);
+               CommandParameter pmismatch("mismatch", "Number", "", "-1.0", "", "", "","",false,false); parameters.push_back(pmismatch);
+               CommandParameter pgapopen("gapopen", "Number", "", "-2.0", "", "", "","",false,false); parameters.push_back(pgapopen);
+               CommandParameter pgapextend("gapextend", "Number", "", "-1.0", "", "", "","",false,false); parameters.push_back(pgapextend);
+               CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
+               CommandParameter poutput("output", "Multiple", "column-lt-square-phylip", "column", "", "", "","phylip-column",false,false,true); parameters.push_back(poutput);
+               CommandParameter pcalc("calc", "Multiple", "nogaps-eachgap-onegap", "onegap", "", "", "","",false,false); parameters.push_back(pcalc);
+               CommandParameter pcountends("countends", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(pcountends);
+               CommandParameter pcompress("compress", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(pcompress);
+               CommandParameter pcutoff("cutoff", "Number", "", "1.0", "", "", "","",false,false,true); parameters.push_back(pcutoff);
+               CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
+               CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
                
                vector<string> myArray;
                for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
@@ -77,7 +64,22 @@ string PairwiseSeqsCommand::getHelpString(){
                exit(1);
        }
 }
-
+//**********************************************************************************************************************
+string PairwiseSeqsCommand::getOutputPattern(string type) {
+    try {
+        string pattern = "";
+        
+        if (type == "phylip") {  pattern = "[filename],[outputtag],dist"; } 
+        else if (type == "column") { pattern = "[filename],dist"; }
+        else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
+        
+        return pattern;
+    }
+    catch(exception& e) {
+        m->errorOut(e, "PairwiseSeqsCommand", "getOutputPattern");
+        exit(1);
+    }
+}
 //**********************************************************************************************************************
 PairwiseSeqsCommand::PairwiseSeqsCommand(){    
        try {
@@ -99,6 +101,7 @@ PairwiseSeqsCommand::PairwiseSeqsCommand(string option)  {
        
                //allow user to run help
                if(option == "help") { help(); abort = true; calledHelp = true; }
+               else if(option == "citation") { citation(); abort = true; calledHelp = true;}
                
                else {
                        vector<string> myArray = setParameters();
@@ -197,6 +200,8 @@ PairwiseSeqsCommand::PairwiseSeqsCommand(string option)  {
                                                        //erase from file list
                                                        fastaFileNames.erase(fastaFileNames.begin()+i);
                                                        i--;
+                                               }else {
+                                                       m->setFastaFile(fastaFileNames[i]);
                                                }
                                        }
                                }
@@ -209,23 +214,26 @@ PairwiseSeqsCommand::PairwiseSeqsCommand(string option)  {
                        // ...at some point should added some additional type checking...
                        string temp;
                        temp = validParameter.validFile(parameters, "match", false);            if (temp == "not found"){       temp = "1.0";                   }
-                       convert(temp, match);  
+                       m->mothurConvert(temp, match);  
                        
                        temp = validParameter.validFile(parameters, "mismatch", false);         if (temp == "not found"){       temp = "-1.0";                  }
-                       convert(temp, misMatch);  
+                       m->mothurConvert(temp, misMatch);  
+            if (misMatch > 0) { m->mothurOut("[ERROR]: mismatch must be negative.\n"); abort=true; }
                        
                        temp = validParameter.validFile(parameters, "gapopen", false);          if (temp == "not found"){       temp = "-2.0";                  }
-                       convert(temp, gapOpen);  
+                       m->mothurConvert(temp, gapOpen);  
+            if (gapOpen > 0) { m->mothurOut("[ERROR]: gapopen must be negative.\n"); abort=true; }
                        
                        temp = validParameter.validFile(parameters, "gapextend", false);        if (temp == "not found"){       temp = "-1.0";                  }
-                       convert(temp, gapExtend); 
+                       m->mothurConvert(temp, gapExtend); 
+            if (gapExtend > 0) { m->mothurOut("[ERROR]: gapextend must be negative.\n"); abort=true; }
                        
                        temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = m->getProcessors();      }
                        m->setProcessors(temp);
-                       convert(temp, processors);
+                       m->mothurConvert(temp, processors);
                        
                        temp = validParameter.validFile(parameters, "cutoff", false);           if(temp == "not found"){        temp = "1.0"; }
-                       convert(temp, cutoff); 
+                       m->mothurConvert(temp, cutoff); 
                        
                        temp = validParameter.validFile(parameters, "countends", false);        if(temp == "not found"){        temp = "T";     }
                        countends = m->isTrue(temp); 
@@ -236,6 +244,7 @@ PairwiseSeqsCommand::PairwiseSeqsCommand(string option)  {
                        align = validParameter.validFile(parameters, "align", false);           if (align == "not found"){      align = "needleman";    }
                        
                        output = validParameter.validFile(parameters, "output", false);         if(output == "not found"){      output = "column"; }
+            if (output=="phylip") { output = "lt"; }
                        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"; }
                        
                        calc = validParameter.validFile(parameters, "calc", false);                     
@@ -244,25 +253,6 @@ PairwiseSeqsCommand::PairwiseSeqsCommand(string option)  {
                                 if (calc == "default")  {  calc = "onegap";  }
                        }
                        m->splitAtDash(calc, Estimators);
-                       
-                       ValidCalculators validCalculator;
-                       if (countends) {
-                               for (int i=0; i<Estimators.size(); i++) {
-                                       if (validCalculator.isValidCalculator("distance", Estimators[i]) == true) { 
-                                               if (Estimators[i] == "nogaps")                  {       distCalculator = new ignoreGaps();      }
-                                               else if (Estimators[i] == "eachgap")    {       distCalculator = new eachGapDist();     }
-                                               else if (Estimators[i] == "onegap")             {       distCalculator = new oneGapDist();      }
-                                       }
-                               }
-                       }else {
-                               for (int i=0; i<Estimators.size(); i++) {
-                                       if (validCalculator.isValidCalculator("distance", Estimators[i]) == true) { 
-                                               if (Estimators[i] == "nogaps")          {       distCalculator = new ignoreGaps();                                      }
-                                               else if (Estimators[i] == "eachgap"){   distCalculator = new eachGapIgnoreTermGapDist();        }
-                                               else if (Estimators[i] == "onegap")     {       distCalculator = new oneGapIgnoreTermGapDist();         }
-                                       }
-                               }
-                       }
                }
                
        }
@@ -277,17 +267,7 @@ int PairwiseSeqsCommand::execute(){
        try {
                if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
                
-               int 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.
-               
-               if(align == "gotoh")                    {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase);                 }
-               else if(align == "needleman")   {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);                                }
-               else if(align == "blast")               {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
-               else if(align == "noalign")             {       alignment = new NoAlign();                                                                                                      }
-               else {
-                       m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
-                       m->mothurOutEndLine();
-                       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
-               }
+               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.
                
                cutoff += 0.005;
                
@@ -306,17 +286,21 @@ int PairwiseSeqsCommand::execute(){
                        int numSeqs = alignDB.getNumSeqs();
                        int startTime = time(NULL);
                        string outputFile = "";
-                               
+                       
+            map<string, string> variables; 
+            variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s]));
                        if (output == "lt") { //does the user want lower triangle phylip formatted file 
-                               outputFile = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + "phylip.dist";
-                               remove(outputFile.c_str()); outputTypes["phylip"].push_back(outputFile);
+                               variables["[outputtag]"] = "phylip";
+                outputFile = getOutputFileName("phylip", variables);
+                               m->mothurRemove(outputFile); outputTypes["phylip"].push_back(outputFile);
                        }else if (output == "column") { //user wants column format
-                               outputFile = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + "dist";
+                               outputFile = getOutputFileName("column", variables);
                                outputTypes["column"].push_back(outputFile);
-                               remove(outputFile.c_str());
+                               m->mothurRemove(outputFile);
                        }else { //assume square
-                               outputFile = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + "square.dist";
-                               remove(outputFile.c_str());
+                variables["[outputtag]"] = "square";
+                outputFile = getOutputFileName("phylip", variables);
+                               m->mothurRemove(outputFile);
                                outputTypes["phylip"].push_back(outputFile);
                        }
                        
@@ -354,11 +338,11 @@ int PairwiseSeqsCommand::execute(){
                                        
                                        driverMPI(start, end, outMPI, cutoff); 
                                        
-                                       if (m->control_pressed) { outputTypes.clear(); MPI_File_close(&outMPI);  remove(outputFile.c_str()); delete distCalculator;  return 0; }
+                                       if (m->control_pressed) { outputTypes.clear(); MPI_File_close(&outMPI);  m->mothurRemove(outputFile); return 0; }
                                
                                        //wait on chidren
                                        for(int i = 1; i < processors; i++) { 
-                                               if (m->control_pressed) { outputTypes.clear();  MPI_File_close(&outMPI);   remove(outputFile.c_str()); delete distCalculator;  return 0; }
+                                               if (m->control_pressed) { outputTypes.clear();  MPI_File_close(&outMPI);   m->mothurRemove(outputFile);  return 0; }
                                                
                                                char buf[5];
                                                MPI_Recv(buf, 5, MPI_CHAR, i, tag, MPI_COMM_WORLD, &status); 
@@ -367,7 +351,7 @@ int PairwiseSeqsCommand::execute(){
                                        //do your part
                                        driverMPI(start, end, outMPI, cutoff); 
                                        
-                                       if (m->control_pressed) { outputTypes.clear();  MPI_File_close(&outMPI);  remove(outputFile.c_str()); delete distCalculator;  return 0; }
+                                       if (m->control_pressed) { outputTypes.clear();  MPI_File_close(&outMPI);  m->mothurRemove(outputFile);  return 0; }
                                
                                        char buf[5];
                                        strcpy(buf, "done"); 
@@ -382,12 +366,12 @@ int PairwiseSeqsCommand::execute(){
                                
                                        //do your part
                                        string outputMyPart;
-                                       unsigned long int mySize;
+                                       unsigned long long mySize;
                                        
                                        if (output != "square"){ driverMPI(start, end, outputFile, mySize); }
                                        else { driverMPI(start, end, outputFile, mySize, output); }
                
-                                       if (m->control_pressed) {  outputTypes.clear();   remove(outputFile.c_str()); delete distCalculator;  return 0; }
+                                       if (m->control_pressed) {  outputTypes.clear();   m->mothurRemove(outputFile);   return 0; }
                                        
                                        int amode=MPI_MODE_APPEND|MPI_MODE_WRONLY|MPI_MODE_CREATE; //
                                        MPI_File outMPI;
@@ -400,9 +384,9 @@ int PairwiseSeqsCommand::execute(){
 
                                        //wait on chidren
                                        for(int b = 1; b < processors; b++) { 
-                                               unsigned long int fileSize;
+                                               unsigned long long fileSize;
                                                
-                                               if (m->control_pressed) { outputTypes.clear();  MPI_File_close(&outMPI);  remove(outputFile.c_str());  delete distCalculator;  return 0; }
+                                               if (m->control_pressed) { outputTypes.clear();  MPI_File_close(&outMPI);  m->mothurRemove(outputFile);   return 0; }
                                                
                                                MPI_Recv(&fileSize, 1, MPI_LONG, b, tag, MPI_COMM_WORLD, &status); 
                                                
@@ -428,11 +412,11 @@ int PairwiseSeqsCommand::execute(){
                                        MPI_File_close(&outMPI);
                                }else { //you are a child process
                                        //do your part
-                                       unsigned long int size;
+                                       unsigned long long size;
                                        if (output != "square"){ driverMPI(start, end, (outputFile + toString(pid) + ".temp"), size); }
                                        else { driverMPI(start, end, (outputFile + toString(pid) + ".temp"), size, output); }
                                        
-                                       if (m->control_pressed) { delete distCalculator;  return 0; }
+                                       if (m->control_pressed) {  return 0; }
                                
                                        //tell parent you are done.
                                        MPI_Send(&size, 1, MPI_LONG, 0, tag, MPI_COMM_WORLD);
@@ -441,7 +425,7 @@ int PairwiseSeqsCommand::execute(){
                        MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
        #else           
                                        
-               #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+               //#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
                        //if you don't need to fork anything
                        if(processors == 1){
                                if (output != "square") {  driver(0, numSeqs, outputFile, cutoff); }
@@ -462,14 +446,14 @@ int PairwiseSeqsCommand::execute(){
                                
                                createProcesses(outputFile); 
                        }
-               #else
+               //#else
                        //ifstream inFASTA;
-                       if (output != "square") {  driver(0, numSeqs, outputFile, cutoff); }
-                       else { driver(0, numSeqs, outputFile, "square");  }
-               #endif
+                       //if (output != "square") {  driver(0, numSeqs, outputFile, cutoff); }
+                       //else { driver(0, numSeqs, outputFile, "square");  }
+               //#endif
                
        #endif
-                       if (m->control_pressed) { outputTypes.clear();  delete distCalculator; remove(outputFile.c_str()); return 0; }
+                       if (m->control_pressed) { outputTypes.clear();   m->mothurRemove(outputFile); return 0; }
                        
                        #ifdef USE_MPI
                                MPI_Comm_rank(MPI_COMM_WORLD, &pid); 
@@ -497,10 +481,8 @@ int PairwiseSeqsCommand::execute(){
                        
                        m->mothurOut("It took " + toString(time(NULL) - startTime) + " to calculate the distances for " + toString(numSeqs) + " sequences."); m->mothurOutEndLine();
                        
-                       if (m->control_pressed) { outputTypes.clear();  delete distCalculator; remove(outputFile.c_str()); return 0; }
+                       if (m->control_pressed) { outputTypes.clear(); m->mothurRemove(outputFile); return 0; }
                }
-                       
-               delete distCalculator;
                
                //set phylip file as new current phylipfile
                string current = "";
@@ -516,7 +498,7 @@ int PairwiseSeqsCommand::execute(){
                }
                
                m->mothurOutEndLine();
-               m->mothurOut("Output File Name: "); m->mothurOutEndLine();
+               m->mothurOut("Output File Names: "); m->mothurOutEndLine();
                for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
                m->mothurOutEndLine();
                
@@ -532,9 +514,11 @@ int PairwiseSeqsCommand::execute(){
 /**************************************************************************************************/
 void PairwiseSeqsCommand::createProcesses(string filename) {
        try {
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
-               int process = 1;
+        int process = 1;
                processIDS.clear();
+        
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+               
                
                //loop through and create all the processes you want
                while (process != processors) {
@@ -564,13 +548,54 @@ void PairwiseSeqsCommand::createProcesses(string filename) {
                        int temp = processIDS[i];
                        wait(&temp);
                }
+#else     
+        //////////////////////////////////////////////////////////////////////////////////////////////////////
+               //Windows version shared memory, so be careful when passing variables through the distanceData struct. 
+               //Above fork() will clone, so memory is separate, but that's not the case with windows, 
+               //that's why the distance calculator was moved inside of the driver to make separate copies.
+               //////////////////////////////////////////////////////////////////////////////////////////////////////
+               
+               vector<pairwiseData*> pDataArray; //[processors-1];
+               DWORD   dwThreadIdArray[processors-1];
+               HANDLE  hThreadArray[processors-1]; 
+               
+               //Create processor-1 worker threads.
+               for( int i=0; i<processors-1; i++ ){
+                       string extension = toString(i) + ".temp";
+
+                       // Allocate memory for thread data.
+                       pairwiseData* tempDist = new pairwiseData((filename+extension), align, "square", Estimators[0], countends, output, alignDB, m, lines[i+1].start, lines[i+1].end, match, misMatch, gapOpen, gapExtend, longestBase, cutoff, i);
+                       pDataArray.push_back(tempDist);
+                       processIDS.push_back(i);
+                       
+                       if (output != "square") { hThreadArray[i] = CreateThread(NULL, 0, MyPairwiseThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);  } 
+            else { hThreadArray[i] = CreateThread(NULL, 0, MyPairwiseSquareThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);  }
+               }
+               
+               //do your part
+               if (output != "square") {  driver(lines[0].start, lines[0].end, filename, cutoff); }
+               else { driver(lines[0].start, lines[0].end, filename, "square"); }
                
-               //append and remove temp files
+               //Wait until all threads have terminated.
+               WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
+               
+               //Close all thread handles and free memory allocations.
+               for(int i=0; i < pDataArray.size(); i++){
+            if (pDataArray[i]->count != (pDataArray[i]->end-pDataArray[i]->start)) {
+                m->mothurOut("[ERROR]: process " + toString(i) + " only processed " + toString(pDataArray[i]->count) + " of " + toString(pDataArray[i]->end-pDataArray[i]->start) + " sequences assigned to it, quitting. \n"); m->control_pressed = true; 
+            }
+                       CloseHandle(hThreadArray[i]);
+                       delete pDataArray[i];
+               }
+
+#endif
+        
+        //append and remove temp files
                for (int i=0;i<processIDS.size();i++) { 
                        m->appendFiles((filename + toString(processIDS[i]) + ".temp"), filename);
-                       remove((filename + toString(processIDS[i]) + ".temp").c_str());
+                       m->mothurRemove((filename + toString(processIDS[i]) + ".temp"));
                }
-#endif
+        
        }
        catch(exception& e) {
                m->errorOut(e, "PairwiseSeqsCommand", "createProcesses");
@@ -584,6 +609,33 @@ int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, fl
        try {
 
                int startTime = time(NULL);
+        
+        Alignment* alignment;
+        if(align == "gotoh")                   {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase);                 }
+               else if(align == "needleman")   {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);                                }
+               else if(align == "blast")               {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
+               else if(align == "noalign")             {       alignment = new NoAlign();                                                                                                      }
+               else {
+                       m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
+                       m->mothurOutEndLine();
+                       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
+               }
+        
+        ValidCalculators validCalculator;
+        Dist* distCalculator;
+        if (countends) {
+            if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) { 
+                if (Estimators[0] == "nogaps")                 {       distCalculator = new ignoreGaps();      }
+                else if (Estimators[0] == "eachgap")   {       distCalculator = new eachGapDist();     }
+                else if (Estimators[0] == "onegap")            {       distCalculator = new oneGapDist();      }
+            }
+        }else {
+            if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) { 
+                if (Estimators[0] == "nogaps")         {       distCalculator = new ignoreGaps();                                      }
+                else if (Estimators[0] == "eachgap"){  distCalculator = new eachGapIgnoreTermGapDist();        }
+                else if (Estimators[0] == "onegap")    {       distCalculator = new oneGapIgnoreTermGapDist();         }
+            }
+        }
                
                //column file
                ofstream outFile(dFileName.c_str(), ios::trunc);
@@ -603,7 +655,7 @@ int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, fl
                        
                        for(int j=0;j<i;j++){
                                
-                               if (m->control_pressed) { outFile.close(); return 0;  }
+                               if (m->control_pressed) { outFile.close(); delete alignment; delete distCalculator; return 0;  }
                                
                                if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
                                        alignment->resize(alignDB.get(i).getUnaligned().length()+1);
@@ -620,10 +672,14 @@ int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, fl
                                seqI.setAligned(alignment->getSeqAAln());
                                seqJ.setAligned(alignment->getSeqBAln());
 
-                               
+                               //cout << seqI.getName() << '\t' << seqJ.getName() << endl;
+                //cout << alignment->getSeqAAln() << endl << alignment->getSeqBAln() << endl;
+                
                                distCalculator->calcDist(seqI, seqJ);
                                double dist = distCalculator->getDist();
-                               
+                
+                //cout << "dist = " << dist << endl;
+                                               
                                if(dist <= cutoff){
                                        if (output == "column") { outFile << alignDB.get(i).getName() << ' ' << alignDB.get(j).getName() << ' ' << dist << endl; }
                                }
@@ -640,6 +696,8 @@ int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, fl
                m->mothurOut(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
                
                outFile.close();
+        delete alignment;
+        delete distCalculator;
                
                return 1;
        }
@@ -654,7 +712,34 @@ int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, st
        try {
 
                int startTime = time(NULL);
+        
+        Alignment* alignment;
+        if(align == "gotoh")                   {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase);                 }
+               else if(align == "needleman")   {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);                                }
+               else if(align == "blast")               {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
+               else if(align == "noalign")             {       alignment = new NoAlign();                                                                                                      }
+               else {
+                       m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
+                       m->mothurOutEndLine();
+                       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
+               }
                
+        ValidCalculators validCalculator;
+        Dist* distCalculator;
+        if (countends) {
+            if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) { 
+                if (Estimators[0] == "nogaps")                 {       distCalculator = new ignoreGaps();      }
+                else if (Estimators[0] == "eachgap")   {       distCalculator = new eachGapDist();     }
+                else if (Estimators[0] == "onegap")            {       distCalculator = new oneGapDist();      }
+            }
+        }else {
+            if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) { 
+                if (Estimators[0] == "nogaps")         {       distCalculator = new ignoreGaps();                                      }
+                else if (Estimators[0] == "eachgap"){  distCalculator = new eachGapIgnoreTermGapDist();        }
+                else if (Estimators[0] == "onegap")    {       distCalculator = new oneGapIgnoreTermGapDist();         }
+            }
+        }
+
                //column file
                ofstream outFile(dFileName.c_str(), ios::trunc);
                outFile.setf(ios::fixed, ios::showpoint);
@@ -672,7 +757,7 @@ int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, st
                        
                        for(int j=0;j<alignDB.getNumSeqs();j++){
                                
-                               if (m->control_pressed) { outFile.close(); return 0;  }
+                               if (m->control_pressed) { outFile.close(); delete alignment; delete distCalculator; return 0;  }
                                
                                if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
                                        alignment->resize(alignDB.get(i).getUnaligned().length()+1);
@@ -705,6 +790,8 @@ int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, st
                m->mothurOut(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
                
                outFile.close();
+        delete alignment;
+        delete distCalculator;
                
                return 1;
        }
@@ -720,14 +807,41 @@ int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, MPI_File& outMPI,
        try {
                MPI_Status status;
                int startTime = time(NULL);
+        
+        Alignment* alignment;
+        if(align == "gotoh")                   {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase);                 }
+               else if(align == "needleman")   {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);                                }
+               else if(align == "blast")               {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
+               else if(align == "noalign")             {       alignment = new NoAlign();                                                                                                      }
+               else {
+                       m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
+                       m->mothurOutEndLine();
+                       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
+               }
                
+        ValidCalculators validCalculator;
+        Dist* distCalculator;
+        if (countends) {
+            if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) { 
+                if (Estimators[0] == "nogaps")                 {       distCalculator = new ignoreGaps();      }
+                else if (Estimators[0] == "eachgap")   {       distCalculator = new eachGapDist();     }
+                else if (Estimators[0] == "onegap")            {       distCalculator = new oneGapDist();      }
+            }
+        }else {
+            if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) { 
+                if (Estimators[0] == "nogaps")         {       distCalculator = new ignoreGaps();                                      }
+                else if (Estimators[0] == "eachgap"){  distCalculator = new eachGapIgnoreTermGapDist();        }
+                else if (Estimators[0] == "onegap")    {       distCalculator = new oneGapIgnoreTermGapDist();         }
+            }
+        }
+
                string outputString = "";
                
                for(int i=startLine;i<endLine;i++){
        
                        for(int j=0;j<i;j++){
                                
-                               if (m->control_pressed) {  return 0;  }
+                               if (m->control_pressed) { delete alignment; delete distCalculator; return 0;  }
                                
                                if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
                                        alignment->resize(alignDB.get(i).getUnaligned().length()+1);
@@ -769,7 +883,8 @@ int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, MPI_File& outMPI,
                        delete buf;
                        
                }
-               
+               delete alignment;
+        delete distCalculator;
                return 1;
        }
        catch(exception& e) {
@@ -779,7 +894,7 @@ int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, MPI_File& outMPI,
 }
 /**************************************************************************************************/
 /////// need to fix to work with calcs and sequencedb
-int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsigned long int& size){
+int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsigned long long& size){
        try {
                MPI_Status status;
                
@@ -791,7 +906,33 @@ int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsi
 
                MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
 
-               
+               Alignment* alignment;
+        if(align == "gotoh")                   {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase);                 }
+               else if(align == "needleman")   {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);                                }
+               else if(align == "blast")               {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
+               else if(align == "noalign")             {       alignment = new NoAlign();                                                                                                      }
+               else {
+                       m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
+                       m->mothurOutEndLine();
+                       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
+               }
+        
+        ValidCalculators validCalculator;
+        Dist* distCalculator;
+        if (countends) {
+            if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) { 
+                if (Estimators[0] == "nogaps")                 {       distCalculator = new ignoreGaps();      }
+                else if (Estimators[0] == "eachgap")   {       distCalculator = new eachGapDist();     }
+                else if (Estimators[0] == "onegap")            {       distCalculator = new oneGapDist();      }
+            }
+        }else {
+            if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) { 
+                if (Estimators[0] == "nogaps")         {       distCalculator = new ignoreGaps();                                      }
+                else if (Estimators[0] == "eachgap"){  distCalculator = new eachGapIgnoreTermGapDist();        }
+                else if (Estimators[0] == "onegap")    {       distCalculator = new oneGapIgnoreTermGapDist();         }
+            }
+        }
+
                
                string outputString = "";
                size = 0;
@@ -808,7 +949,7 @@ int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsi
                        
                        for(int j=0;j<i;j++){
                                
-                               if (m->control_pressed) {  return 0;  }
+                               if (m->control_pressed) { delete alignment; delete distCalculator; return 0;  }
                                
                                if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
                                        alignment->resize(alignDB.get(i).getUnaligned().length()+1);
@@ -845,6 +986,8 @@ int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsi
                }
                
                MPI_File_close(&outMPI);
+        delete alignment;
+        delete distCalculator;
                
                return 1;
        }
@@ -855,7 +998,7 @@ int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsi
 }
 /**************************************************************************************************/
 /////// need to fix to work with calcs and sequencedb
-int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsigned long int& size, string square){
+int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsigned long long& size, string square){
        try {
                MPI_Status status;
                
@@ -867,8 +1010,33 @@ int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsi
 
                MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
                
+               Alignment* alignment;
+        if(align == "gotoh")                   {       alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase);                 }
+               else if(align == "needleman")   {       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);                                }
+               else if(align == "blast")               {       alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);            }
+               else if(align == "noalign")             {       alignment = new NoAlign();                                                                                                      }
+               else {
+                       m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
+                       m->mothurOutEndLine();
+                       alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
+               }
                
-               
+        ValidCalculators validCalculator;
+        Dist* distCalculator;
+        if (countends) {
+            if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) { 
+                if (Estimators[0] == "nogaps")                 {       distCalculator = new ignoreGaps();      }
+                else if (Estimators[0] == "eachgap")   {       distCalculator = new eachGapDist();     }
+                else if (Estimators[0] == "onegap")            {       distCalculator = new oneGapDist();      }
+            }
+        }else {
+            if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) { 
+                if (Estimators[0] == "nogaps")         {       distCalculator = new ignoreGaps();                                      }
+                else if (Estimators[0] == "eachgap"){  distCalculator = new eachGapIgnoreTermGapDist();        }
+                else if (Estimators[0] == "onegap")    {       distCalculator = new oneGapIgnoreTermGapDist();         }
+            }
+        }
+        
                string outputString = "";
                size = 0;
                
@@ -884,7 +1052,7 @@ int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsi
                        
                        for(int j=0;j<alignDB.getNumSeqs();j++){
                                
-                               if (m->control_pressed) {  return 0;  }
+                               if (m->control_pressed) {  delete alignment; return 0;  }
                                
                                if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
                                        alignment->resize(alignDB.get(i).getUnaligned().length()+1);
@@ -922,6 +1090,7 @@ int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsi
                
                MPI_File_close(&outMPI);
                
+        delete alignment;
                return 1;
        }
        catch(exception& e) {