]> git.donarmstrong.com Git - mothur.git/blobdiff - aligncommand.cpp
fixed get.lienage and remove.lineage bug with confidence scores that are floats....
[mothur.git] / aligncommand.cpp
index b4af28a2990a44d69527d4c5cb7e5984d8b73033..a820542e42f4d7a7b9956b194f790555e634536a 100644 (file)
@@ -64,7 +64,7 @@ vector<string> AlignCommand::getRequiredFiles(){
 //**********************************************************************************************************************
 AlignCommand::AlignCommand(){  
        try {
-               //initialize outputTypes
+               abort = true; calledHelp = true; 
                vector<string> tempOutNames;
                outputTypes["fasta"] = tempOutNames;
                outputTypes["alignreport"] = tempOutNames;
@@ -78,10 +78,10 @@ AlignCommand::AlignCommand(){
 //**********************************************************************************************************************
 AlignCommand::AlignCommand(string option)  {
        try {
-               abort = false;
+               abort = false; calledHelp = false;   
        
                //allow user to run help
-               if(option == "help") { help(); abort = true; }
+               if(option == "help") { help(); abort = true; calledHelp = true;}
                
                else {
                        
@@ -137,8 +137,8 @@ AlignCommand::AlignCommand(string option)  {
                        }else if (templateFileName == "not open") { abort = true; }     
                        
                        candidateFileName = validParameter.validFile(parameters, "candidate", false);
-                       if (candidateFileName == "not found") { m->mothurOut("candidate is a required parameter for the align.seqs command."); m->mothurOutEndLine(); abort = true;  }
-                       else { 
+                       if (candidateFileName == "not found") { 
+                       }else { 
                                m->splitAtDash(candidateFileName, candidateFileNames);
                                
                                //go through files and make sure they are good, if not, then disregard them
@@ -153,30 +153,34 @@ AlignCommand::AlignCommand(string option)  {
        
                                        int ableToOpen;
                                        ifstream in;
-
                                        ableToOpen = m->openInputFile(candidateFileNames[i], in, "noerror");
-                               
+                                       in.close();     
+                                       
                                        //if you can't open it, try default location
                                        if (ableToOpen == 1) {
                                                if (m->getDefaultPath() != "") { //default path is set
                                                        string tryPath = m->getDefaultPath() + m->getSimpleName(candidateFileNames[i]);
                                                        m->mothurOut("Unable to open " + candidateFileNames[i] + ". Trying default " + tryPath); m->mothurOutEndLine();
-                                                       ableToOpen = m->openInputFile(tryPath, in, "noerror");
+                                                       ifstream in2;
+                                                       ableToOpen = m->openInputFile(tryPath, in2, "noerror");
+                                                       in2.close();
                                                        candidateFileNames[i] = tryPath;
                                                }
                                        }
                                        
-                                       //if you can't open it, try default location
+                                       //if you can't open it, try output location
                                        if (ableToOpen == 1) {
                                                if (m->getOutputDir() != "") { //default path is set
                                                        string tryPath = m->getOutputDir() + m->getSimpleName(candidateFileNames[i]);
                                                        m->mothurOut("Unable to open " + candidateFileNames[i] + ". Trying output directory " + tryPath); m->mothurOutEndLine();
-                                                       ableToOpen = m->openInputFile(tryPath, in, "noerror");
+                                                       ifstream in2;
+                                                       ableToOpen = m->openInputFile(tryPath, in2, "noerror");
+                                                       in2.close();
                                                        candidateFileNames[i] = tryPath;
                                                }
                                        }
                                        
-                                       in.close();                                     
+                                                                       
 
                                        if (ableToOpen == 1) { 
                                                m->mothurOut("Unable to open " + candidateFileNames[i] + ". It will be disregarded."); m->mothurOutEndLine(); 
@@ -274,7 +278,7 @@ void AlignCommand::help(){
 
 int AlignCommand::execute(){
        try {
-               if (abort == true) {    return 0;       }
+               if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 
                templateDB = new AlignmentDB(templateFileName, search, kmerSize, gapOpen, gapExtend, match, misMatch);
                int longestBase = templateDB->getLongestBase();
@@ -305,7 +309,7 @@ int AlignCommand::execute(){
                        int start = time(NULL);
                
 #ifdef USE_MPI 
-                               int pid, end, numSeqsPerProcessor; 
+                               int pid, numSeqsPerProcessor; 
                                int tag = 2001;
                                vector<unsigned long int> MPIPos;
                                MPIWroteAccnos = false;
@@ -418,62 +422,12 @@ int AlignCommand::execute(){
        #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
                        if(processors == 1){
                                numFastaSeqs = driver(lines[0], alignFileName, reportFileName, accnosFileName, candidateFileNames[s]);
-                               
-                               if (m->control_pressed) { remove(accnosFileName.c_str()); remove(alignFileName.c_str()); remove(reportFileName.c_str()); outputTypes.clear(); return 0; }
-                               
-                               //delete accnos file if its blank else report to user
-                               if (m->isBlank(accnosFileName)) {  remove(accnosFileName.c_str());  hasAccnos = false; }
-                               else { 
-                                       m->mothurOut("Some of you sequences generated alignments that eliminated too many bases, a list is provided in " + accnosFileName + ".");
-                                       if (!flip) {
-                                               m->mothurOut(" If you set the flip parameter to true mothur will try aligning the reverse compliment as well."); 
-                                       }else{  m->mothurOut(" If the reverse compliment proved to be better it was reported.");  }
-                                       m->mothurOutEndLine();
-                               }
                        }else{
-                               processIDS.resize(0);
-                               
                                numFastaSeqs = createProcesses(alignFileName, reportFileName, accnosFileName, candidateFileNames[s]); 
-                               
-                               rename((alignFileName + toString(processIDS[0]) + ".temp").c_str(), alignFileName.c_str());
-                               rename((reportFileName + toString(processIDS[0]) + ".temp").c_str(), reportFileName.c_str());
-                               
-                               //append alignment and report files
-                               for(int i=1;i<processors;i++){
-                                       appendAlignFiles((alignFileName + toString(processIDS[i]) + ".temp"), alignFileName);
-                                       remove((alignFileName + toString(processIDS[i]) + ".temp").c_str());
-                                       
-                                       appendReportFiles((reportFileName + toString(processIDS[i]) + ".temp"), reportFileName);
-                                       remove((reportFileName + toString(processIDS[i]) + ".temp").c_str());
-                               }
-                               
-                               vector<string> nonBlankAccnosFiles;
-                               //delete blank accnos files generated with multiple processes
-                               for(int i=0;i<processors;i++){  
-                                       if (!(m->isBlank(accnosFileName + toString(processIDS[i]) + ".temp"))) {
-                                               nonBlankAccnosFiles.push_back(accnosFileName + toString(processIDS[i]) + ".temp");
-                                       }else { remove((accnosFileName + toString(processIDS[i]) + ".temp").c_str());  }
-                               }
-                               
-                               //append accnos files
-                               if (nonBlankAccnosFiles.size() != 0) { 
-                                       rename(nonBlankAccnosFiles[0].c_str(), accnosFileName.c_str());
-                                       
-                                       for (int h=1; h < nonBlankAccnosFiles.size(); h++) {
-                                               appendAlignFiles(nonBlankAccnosFiles[h], accnosFileName);
-                                               remove(nonBlankAccnosFiles[h].c_str());
-                                       }
-                                       m->mothurOut("Some of you sequences generated alignments that eliminated too many bases, a list is provided in " + accnosFileName + ".");
-                                       if (!flip) {
-                                               m->mothurOut(" If you set the flip parameter to true mothur will try aligning the reverse compliment as well."); 
-                                       }else{  m->mothurOut(" If the reverse compliment proved to be better it was reported.");  }
-                                       m->mothurOutEndLine();
-                               }else{ hasAccnos = false;  }
-                               
-                               if (m->control_pressed) { remove(accnosFileName.c_str()); remove(alignFileName.c_str()); remove(reportFileName.c_str()); outputTypes.clear(); return 0; }
                        }
        #else
                        numFastaSeqs = driver(lines[0], alignFileName, reportFileName, accnosFileName, candidateFileNames[s]);
+       #endif
                        
                        if (m->control_pressed) { remove(accnosFileName.c_str()); remove(alignFileName.c_str()); remove(reportFileName.c_str()); outputTypes.clear();  return 0; }
                        
@@ -486,8 +440,6 @@ int AlignCommand::execute(){
                                }else{  m->mothurOut(" If the reverse compliment proved to be better it was reported.");  }
                                m->mothurOutEndLine();
                        }
-                       
-       #endif
 
 #endif         
 
@@ -511,6 +463,12 @@ int AlignCommand::execute(){
                        m->mothurOutEndLine();
                }
                
+               //set align file as new current fastafile
+               string currentFasta = "";
+               itTypes = outputTypes.find("fasta");
+               if (itTypes != outputTypes.end()) {
+                       if ((itTypes->second).size() != 0) { currentFasta = (itTypes->second)[0]; m->setFastaFile(currentFasta); }
+               }
                
                m->mothurOutEndLine();
                m->mothurOut("Output File Names: "); m->mothurOutEndLine();
@@ -550,7 +508,8 @@ int AlignCommand::driver(linePair* filePos, string alignFName, string reportFNam
                        if (m->control_pressed) {  return 0; }
                        
                        Sequence* candidateSeq = new Sequence(inFASTA);  m->gobble(inFASTA);
-                       
+                       report.setCandidate(candidateSeq);
+
                        int origNumBases = candidateSeq->getNumBases();
                        string originalUnaligned = candidateSeq->getUnaligned();
                        int numBasesNeeded = origNumBases * threshold;
@@ -562,10 +521,11 @@ int AlignCommand::driver(linePair* filePos, string alignFName, string reportFNam
                                                                
                                Sequence temp = templateDB->findClosestSequence(candidateSeq);
                                Sequence* templateSeq = &temp;
-                               
+                       
                                float searchScore = templateDB->getSearchScore();
                                                                
                                Nast* nast = new Nast(alignment, candidateSeq, templateSeq);
+               
                                Sequence* copy;
                                
                                Nast* nast2;
@@ -580,6 +540,7 @@ int AlignCommand::driver(linePair* filePos, string alignFName, string reportFNam
                                        string wasBetter =  "";
                                        //if the user wants you to try the reverse
                                        if (flip) {
+                               
                                                //get reverse compliment
                                                copy = new Sequence(candidateSeq->getName(), originalUnaligned);
                                                copy->reverseComplement();
@@ -611,7 +572,6 @@ int AlignCommand::driver(linePair* filePos, string alignFName, string reportFNam
                                        accnosFile << candidateSeq->getName() << wasBetter << endl;
                                }
                                
-                               report.setCandidate(candidateSeq);
                                report.setTemplate(templateSeq);
                                report.setSearchParameters(search, searchScore);
                                report.setAlignmentParameters(align, alignment);
@@ -699,7 +659,8 @@ int AlignCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& align
                        istringstream iss (tempBuf,istringstream::in);
 
                        Sequence* candidateSeq = new Sequence(iss);  
-       
+                       report.setCandidate(candidateSeq);
+
                        int origNumBases = candidateSeq->getNumBases();
                        string originalUnaligned = candidateSeq->getUnaligned();
                        int numBasesNeeded = origNumBases * threshold;
@@ -748,6 +709,7 @@ int AlignCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& align
                                                        delete nast;
                                                        nast = nast2;
                                                        needToDeleteCopy = true;
+                                                       wasBetter = "\treverse complement produced a better alignment, so mothur used the reverse complement.";
                                                }else{  
                                                        wasBetter = "\treverse complement did NOT produce a better alignment, please check sequence.";
                                                        delete nast2;
@@ -769,7 +731,6 @@ int AlignCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& align
                                        MPIWroteAccnos = true;
                                }
                                
-                               report.setCandidate(candidateSeq);
                                report.setTemplate(templateSeq);
                                report.setSearchParameters(search, searchScore);
                                report.setAlignmentParameters(align, alignment);
@@ -820,7 +781,8 @@ int AlignCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& align
 int AlignCommand::createProcesses(string alignFileName, string reportFileName, string accnosFName, string filename) {
        try {
 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
-               int process = 0;
+               processIDS.resize(0);
+               int process = 1;
                int num = 0;
                //              processIDS.resize(0);
                
@@ -842,21 +804,57 @@ int AlignCommand::createProcesses(string alignFileName, string reportFileName, s
                                out.close();
                                
                                exit(0);
-                       }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
+                       }else { 
+                               m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); 
+                               for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
+                               exit(0);
+                       }
                }
                
+               //do my part
+               num = driver(lines[0], alignFileName, reportFileName, accnosFName, filename);
+               
                //force parent to wait until all the processes are done
                for (int i=0;i<processors;i++) { 
                        int temp = processIDS[i];
                        wait(&temp);
                }
                
+               vector<string> nonBlankAccnosFiles;
+               if (!(m->isBlank(accnosFName))) { nonBlankAccnosFiles.push_back(accnosFName); }
+               else { remove(accnosFName.c_str()); } //remove so other files can be renamed to it
+                       
                for (int i = 0; i < processIDS.size(); i++) {
                        ifstream in;
                        string tempFile =  alignFileName + toString(processIDS[i]) + ".num.temp";
                        m->openInputFile(tempFile, in);
                        if (!in.eof()) { int tempNum = 0; in >> tempNum; num += tempNum; }
                        in.close(); remove(tempFile.c_str());
+                       
+                       appendAlignFiles((alignFileName + toString(processIDS[i]) + ".temp"), alignFileName);
+                       remove((alignFileName + toString(processIDS[i]) + ".temp").c_str());
+                       
+                       appendReportFiles((reportFileName + toString(processIDS[i]) + ".temp"), reportFileName);
+                       remove((reportFileName + toString(processIDS[i]) + ".temp").c_str());
+                       
+                       if (!(m->isBlank(accnosFName + toString(processIDS[i]) + ".temp"))) {
+                               nonBlankAccnosFiles.push_back(accnosFName + toString(processIDS[i]) + ".temp");
+                       }else { remove((accnosFName + toString(processIDS[i]) + ".temp").c_str());  }
+                       
+               }
+               
+               //append accnos files
+               if (nonBlankAccnosFiles.size() != 0) { 
+                       rename(nonBlankAccnosFiles[0].c_str(), accnosFName.c_str());
+                       
+                       for (int h=1; h < nonBlankAccnosFiles.size(); h++) {
+                               appendAlignFiles(nonBlankAccnosFiles[h], accnosFName);
+                               remove(nonBlankAccnosFiles[h].c_str());
+                       }
+               }else { //recreate the accnosfile if needed
+                       ofstream out;
+                       m->openOutputFile(accnosFName, out);
+                       out.close();
                }
                
                return num;