]> git.donarmstrong.com Git - mothur.git/blobdiff - preclustercommand.cpp
fix trim.seqs / qualscores bug
[mothur.git] / preclustercommand.cpp
index c5b66a81eaa7ec7381af79c4170f1efaf41e91d1..14409f5a0a6aa86a166f42ff6f99423777237759 100644 (file)
 //**********************************************************************************************************************
 inline bool comparePriority(seqPNode first, seqPNode second) {  return (first.numIdentical > second.numIdentical); }
 //**********************************************************************************************************************
+vector<string> PreClusterCommand::getValidParameters(){        
+       try {
+               string Array[] =  {"fasta", "name", "diffs", "outputdir","inputdir"};
+               vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+               return myArray;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "PreClusterCommand", "getValidParameters");
+               exit(1);
+       }
+}
+//**********************************************************************************************************************
+PreClusterCommand::PreClusterCommand(){        
+       try {
+               abort = true; calledHelp = true; 
+               vector<string> tempOutNames;
+               outputTypes["fasta"] = tempOutNames;
+               outputTypes["name"] = tempOutNames;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "PreClusterCommand", "PreClusterCommand");
+               exit(1);
+       }
+}
+//**********************************************************************************************************************
+vector<string> PreClusterCommand::getRequiredParameters(){     
+       try {
+               string Array[] =  {"fasta"};
+               vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+               return myArray;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "PreClusterCommand", "getRequiredParameters");
+               exit(1);
+       }
+}
+//**********************************************************************************************************************
+vector<string> PreClusterCommand::getRequiredFiles(){  
+       try {
+               vector<string> myArray;
+               return myArray;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "PreClusterCommand", "getRequiredFiles");
+               exit(1);
+       }
+}
+//**********************************************************************************************************************
 
 PreClusterCommand::PreClusterCommand(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 {
                        //valid paramters for this command
@@ -36,6 +84,11 @@ PreClusterCommand::PreClusterCommand(string option) {
                                if (validParameter.isValidParameter(it2->first, myArray, it2->second) != true) {  abort = true;  }
                        }
                        
+                       //initialize outputTypes
+                       vector<string> tempOutNames;
+                       outputTypes["fasta"] = tempOutNames;
+                       outputTypes["name"] = tempOutNames;
+               
                        //if the user changes the input directory command factory will send this info to us in the output parameter 
                        string inputDir = validParameter.validFile(parameters, "inputdir", false);              
                        if (inputDir == "not found"){   inputDir = "";          }
@@ -44,7 +97,7 @@ PreClusterCommand::PreClusterCommand(string option) {
                                it = parameters.find("fasta");
                                //user has given a template file
                                if(it != parameters.end()){ 
-                                       path = hasPath(it->second);
+                                       path = m->hasPath(it->second);
                                        //if the user has not given a path then, add inputdir. else leave path alone.
                                        if (path == "") {       parameters["fasta"] = inputDir + it->second;            }
                                }
@@ -52,7 +105,7 @@ PreClusterCommand::PreClusterCommand(string option) {
                                it = parameters.find("name");
                                //user has given a template file
                                if(it != parameters.end()){ 
-                                       path = hasPath(it->second);
+                                       path = m->hasPath(it->second);
                                        //if the user has not given a path then, add inputdir. else leave path alone.
                                        if (path == "") {       parameters["name"] = inputDir + it->second;             }
                                }
@@ -66,7 +119,7 @@ PreClusterCommand::PreClusterCommand(string option) {
                        //if the user changes the output directory command factory will send this info to us in the output parameter 
                        outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  
                                outputDir = ""; 
-                               outputDir += hasPath(fastafile); //if user entered a file with a path then preserve it  
+                               outputDir += m->hasPath(fastafile); //if user entered a file with a path then preserve it       
                        }
 
                        //check for optional parameter and set defaults
@@ -113,7 +166,7 @@ void PreClusterCommand::help(){
 int PreClusterCommand::execute(){
        try {
                
-               if (abort == true) { return 0; }
+               if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
                
                int start = time(NULL);
                
@@ -125,76 +178,72 @@ int PreClusterCommand::execute(){
                if (numSeqs == 0) { m->mothurOut("Error reading fasta file...please correct."); m->mothurOutEndLine(); return 0;  }
                if (diffs > length) { m->mothurOut("Error: diffs is greater than your sequence length."); m->mothurOutEndLine(); return 0;  }
                
-               string fileroot = outputDir + getRootName(getSimpleName(fastafile));
-               string newFastaFile = fileroot + "precluster" + getExtension(fastafile);
-               string newNamesFile = fileroot + "precluster.names";
-               ofstream outFasta;
-               ofstream outNames;
-               
-               openOutputFile(newFastaFile, outFasta);
-               openOutputFile(newNamesFile, outNames);
-
+               //clear sizes since you only needed this info to build the alignSeqs seqPNode structs
+//             sizes.clear();
+       
                //sort seqs by number of identical seqs
-               alignSeqs.sort(comparePriority);
-
+               sort(alignSeqs.begin(), alignSeqs.end(), comparePriority);
+       
                int count = 0;
-               int i = 0;
+
                //think about running through twice...
-               list<seqPNode>::iterator itList;
-               list<seqPNode>::iterator itList2;
-               for (itList = alignSeqs.begin(); itList != alignSeqs.end();) {
+               for (int i = 0; i < numSeqs; i++) {
                        
-                       //try to merge it with all smaller seqs
-                       for (itList2 = alignSeqs.begin(); itList2 != alignSeqs.end();) {
-                               
-                               if (m->control_pressed) { outFasta.close(); outNames.close(); remove(newFastaFile.c_str()); remove(newNamesFile.c_str());  return 0; }
-                               
+                       //are you active
+                       //                      itActive = active.find(alignSeqs[i].seq.getName());
                        
-                               if (itList->seq.getName() != itList2->seq.getName()) { //you don't want to merge with yourself
-                                       //are you within "diff" bases
+                       if (alignSeqs[i].active) {  //this sequence has not been merged yet
+                               
+                               //try to merge it with all smaller seqs
+                               for (int j = i+1; j < numSeqs; j++) {
                                        
-                                       int mismatch = calcMisMatches((*itList).seq.getAligned(), (*itList2).seq.getAligned());
-
-                                       if (mismatch <= diffs) {
-                                               //merge
-                                               (*itList).names += ',' + (*itList2).names;
-                                               (*itList).numIdentical += (*itList2).numIdentical;
+                                       if (m->control_pressed) { return 0; }
+                                       
+                                       if (alignSeqs[j].active) {  //this sequence has not been merged yet
+                                               //are you within "diff" bases
+                                               int mismatch = calcMisMatches(alignSeqs[i].seq.getAligned(), alignSeqs[j].seq.getAligned());
                                                
-                                               itList2 = alignSeqs.erase(itList2); //itList2--;
-                                               count++;
-                                       }else{ itList2++; }
-                               }else{ itList2++; }
-
-                       }
+                                               if (mismatch <= diffs) {
+                                                       //merge
+                                                       alignSeqs[i].names += ',' + alignSeqs[j].names;
+                                                       alignSeqs[i].numIdentical += alignSeqs[j].numIdentical;
 
-                       //ouptut this sequence
-                       printData(outFasta, outNames, (*itList));
-                       
-                       //remove sequence
-                       itList = alignSeqs.erase(itList); 
-                       
-                       i++;
+                                                       alignSeqs[j].active = 0;
+                                                       alignSeqs[j].numIdentical = 0;
+                                                       count++;
+                                               }
+                                       }//end if j active
+                               }//end if i != j
                        
+                       //remove from active list 
+                               alignSeqs[i].active = 0;
+                               
+                       }//end if active i
                        if(i % 100 == 0)        { m->mothurOut(toString(i) + "\t" + toString(numSeqs - count) + "\t" + toString(count)); m->mothurOutEndLine(); }
                }
-                               
-               if(i % 100 != 0)        { m->mothurOut(toString(i) + "\t" + toString(numSeqs - count) + "\t" + toString(count)); m->mothurOutEndLine(); }
                
-               outFasta.close();
-               outNames.close();
+               if(numSeqs % 100 != 0)  { m->mothurOut(toString(numSeqs) + "\t" + toString(numSeqs - count) + "\t" + toString(count)); m->mothurOutEndLine();   }
+       
                
-               if (m->control_pressed) {  remove(newFastaFile.c_str()); remove(newNamesFile.c_str());  return 0; }
-
-               m->mothurOut("It took " + toString(time(NULL) - start) + " secs to cluster " + toString(numSeqs) + " sequences.");
-               m->mothurOut("Total number of sequences before precluster was " + toString(numSeqs) + "."); m->mothurOutEndLine();
-               m->mothurOut("pre.cluster removed " + toString(count) + " sequences."); m->mothurOutEndLine(); 
+               string fileroot = outputDir + m->getRootName(m->getSimpleName(fastafile));
+               
+               string newFastaFile = fileroot + "precluster" + m->getExtension(fastafile);
+               string newNamesFile = fileroot + "precluster.names";
+               
+               if (m->control_pressed) { return 0; }
+               
+               m->mothurOut("Total number of sequences before precluster was " + toString(alignSeqs.size()) + "."); m->mothurOutEndLine();
+               m->mothurOut("pre.cluster removed " + toString(count) + " sequences."); m->mothurOutEndLine(); m->mothurOutEndLine(); 
+               printData(newFastaFile, newNamesFile);
+               
+               m->mothurOut("It took " + toString(time(NULL) - start) + " secs to cluster " + toString(numSeqs) + " sequences."); m->mothurOutEndLine(); 
                
                if (m->control_pressed) { remove(newFastaFile.c_str()); remove(newNamesFile.c_str()); return 0; }
                
                m->mothurOutEndLine();
                m->mothurOut("Output File Names: "); m->mothurOutEndLine();
-               m->mothurOut(newFastaFile); m->mothurOutEndLine();      
-               m->mothurOut(newNamesFile); m->mothurOutEndLine();      
+               m->mothurOut(newFastaFile); m->mothurOutEndLine();      outputNames.push_back(newFastaFile); outputTypes["fasta"].push_back(newFastaFile);
+               m->mothurOut(newNamesFile); m->mothurOutEndLine();      outputNames.push_back(newNamesFile); outputTypes["name"].push_back(newNamesFile);
                m->mothurOutEndLine();
 
                return 0;
@@ -215,8 +264,8 @@ int PreClusterCommand::readFASTA(){
                //ifstream inNames;
                ifstream inFasta;
                
-               //openInputFile(namefile, inNames);
-               openInputFile(fastafile, inFasta);
+               //m->openInputFile(namefile, inNames);
+               m->openInputFile(fastafile, inFasta);
                
                //string firstCol, secondCol, nameString;
                length = 0;
@@ -228,14 +277,14 @@ int PreClusterCommand::readFASTA(){
                        //inNames >> firstCol >> secondCol;
                        //nameString = secondCol;
                        
-                       //gobble(inNames);
+                       //m->gobble(inNames);
                        //int size = 1;
                        //while (secondCol.find_first_of(',') != -1) { 
                        //      size++;
                        //      secondCol = secondCol.substr(secondCol.find_first_of(',')+1, secondCol.length());
                        //}
                        
-                       Sequence seq(inFasta);  gobble(inFasta);
+                       Sequence seq(inFasta);  m->gobble(inFasta);
                        
                        if (seq.getName() != "") {  //can get "" if commented line is at end of fasta file
                                if (namefile != "") {
@@ -287,10 +336,25 @@ int PreClusterCommand::calcMisMatches(string seq1, string seq2){
 
 /**************************************************************************************************/
 
-void PreClusterCommand::printData(ofstream& outFasta, ofstream& outNames, seqPNode thisSeq){
+void PreClusterCommand::printData(string newfasta, string newname){
        try {
-               thisSeq.seq.printSequence(outFasta); 
-               outNames << thisSeq.seq.getName() << '\t' << thisSeq.names << endl;
+               ofstream outFasta;
+               ofstream outNames;
+               
+               m->openOutputFile(newfasta, outFasta);
+               m->openOutputFile(newname, outNames);
+                               
+               
+               for (int i = 0; i < alignSeqs.size(); i++) {
+                       if (alignSeqs[i].numIdentical != 0) {
+                               alignSeqs[i].seq.printSequence(outFasta); 
+                               outNames << alignSeqs[i].seq.getName() << '\t' << alignSeqs[i].names << endl;
+                       }
+               }
+               
+               outFasta.close();
+               outNames.close();
+               
        }
        catch(exception& e) {
                m->errorOut(e, "PreClusterCommand", "printData");
@@ -302,11 +366,11 @@ void PreClusterCommand::printData(ofstream& outFasta, ofstream& outNames, seqPNo
 void PreClusterCommand::readNameFile(){
        try {
                ifstream in;
-               openInputFile(namefile, in);
+               m->openInputFile(namefile, in);
                string firstCol, secondCol;
                                
                while (!in.eof()) {
-                       in >> firstCol >> secondCol; gobble(in);
+                       in >> firstCol >> secondCol; m->gobble(in);
                        names[firstCol] = secondCol;
                        int size = 1;