]> git.donarmstrong.com Git - mothur.git/blobdiff - preclustercommand.cpp
changes while testing
[mothur.git] / preclustercommand.cpp
index a1d5f48bdd77c81bb12d0095628922202a4bcacd..bd3f172c2bcaa3e13ced0b5c57d92e39c9073969 100644 (file)
@@ -17,7 +17,7 @@ vector<string> PreClusterCommand::setParameters(){
                CommandParameter pname("name", "InputTypes", "", "", "NameCount", "none", "none","name",false,false,true); parameters.push_back(pname);
         CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "none", "none","count",false,false,true); parameters.push_back(pcount);
                CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "none", "none","",false,false,true); parameters.push_back(pgroup);
-               CommandParameter pdiffs("diffs", "Number", "", "0", "", "", "","",false,false,true); parameters.push_back(pdiffs);
+               CommandParameter pdiffs("diffs", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pdiffs);
                CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
         CommandParameter ptopdown("topdown", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(ptopdown);
 
@@ -44,7 +44,7 @@ string PreClusterCommand::getHelpString(){
                helpString += "The group parameter allows you to provide a group file so you can cluster by group. \n";
         helpString += "The count parameter allows you to provide a count file so you can cluster by group. \n";
                helpString += "The diffs parameter allows you to specify maximum number of mismatched bases allowed between sequences in a grouping. The default is 1.\n";
-        helpString += "The topdown parameter allows you to specify whether to cluster from largest abundance to smallest or smallest to largest.  Default=T, meanging largest to smallest.\n";
+        helpString += "The topdown parameter allows you to specify whether to cluster from largest abundance to smallest or smallest to largest.  Default=T, meaning largest to smallest.\n";
                helpString += "The pre.cluster command should be in the following format: \n";
                helpString += "pre.cluster(fasta=yourFastaFile, names=yourNamesFile, diffs=yourMaxDiffs) \n";
                helpString += "Example pre.cluster(fasta=amazon.fasta, diffs=2).\n";
@@ -192,7 +192,7 @@ PreClusterCommand::PreClusterCommand(string option) {
                        else if (countfile == "not open") { abort = true; countfile =  ""; }    
                        else {   
                 m->setCountTableFile(countfile); 
-                ct.readTable(countfile);
+                ct.readTable(countfile, true);
                 if (ct.hasGroupInfo()) { bygroup = true; }
                 else { bygroup = false;  }
             }
@@ -560,47 +560,88 @@ int PreClusterCommand::process(string newMapFile){
                int count = 0;
                int numSeqs = alignSeqs.size();
                
-               //think about running through twice...
-               for (int i = 0; i < numSeqs; i++) {
-                       
-                       //are you active
-                       //                      itActive = active.find(alignSeqs[i].seq.getName());
-                       
-                       if (alignSeqs[i].active) {  //this sequence has not been merged yet
-                               
-                               string chunk = alignSeqs[i].seq.getName() + "\t" + toString(alignSeqs[i].numIdentical) + "\t" + toString(0) + "\t" + alignSeqs[i].seq.getAligned() + "\n";
-                               
-                               //try to merge it with all smaller seqs
-                               for (int j = i+1; j < numSeqs; j++) {
-                                       
-                                       if (m->control_pressed) { out.close(); 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());
-                                               
-                                               if (mismatch <= diffs) {
-                                                       //merge
-                                                       alignSeqs[i].names += ',' + alignSeqs[j].names;
-                                                       alignSeqs[i].numIdentical += alignSeqs[j].numIdentical;
-                                                       
-                                                       chunk += alignSeqs[j].seq.getName() + "\t" + toString(alignSeqs[j].numIdentical) + "\t" + toString(mismatch) + "\t" + alignSeqs[j].seq.getAligned() + "\n";
-                                                       
-                                                       alignSeqs[j].active = 0;
-                                                       alignSeqs[j].numIdentical = 0;
-                                                       count++;
-                                               }
-                                       }//end if j active
-                               }//end for loop j
-                               
-                               //remove from active list 
-                               alignSeqs[i].active = 0;
-                               
-                               out << "ideal_seq_" << (i+1) << '\t' << alignSeqs[i].numIdentical << endl << chunk << endl;;
-                               
-                       }//end if active i
-                       if(i % 100 == 0)        { m->mothurOut(toString(i) + "\t" + toString(numSeqs - count) + "\t" + toString(count)); m->mothurOutEndLine(); }
-               }
+        if (topdown) {
+            //think about running through twice...
+            for (int i = 0; i < numSeqs; i++) {
+                
+                if (alignSeqs[i].active) {  //this sequence has not been merged yet
+                    
+                    string chunk = alignSeqs[i].seq.getName() + "\t" + toString(alignSeqs[i].numIdentical) + "\t" + toString(0) + "\t" + alignSeqs[i].seq.getAligned() + "\n";
+                    
+                    //try to merge it with all smaller seqs
+                    for (int j = i+1; j < numSeqs; j++) {
+                        
+                        if (m->control_pressed) { out.close(); 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());
+                            
+                            if (mismatch <= diffs) {
+                                //merge
+                                alignSeqs[i].names += ',' + alignSeqs[j].names;
+                                alignSeqs[i].numIdentical += alignSeqs[j].numIdentical;
+                                
+                                chunk += alignSeqs[j].seq.getName() + "\t" + toString(alignSeqs[j].numIdentical) + "\t" + toString(mismatch) + "\t" + alignSeqs[j].seq.getAligned() + "\n";
+                                
+                                alignSeqs[j].active = 0;
+                                alignSeqs[j].numIdentical = 0;
+                                count++;
+                            }
+                        }//end if j active
+                    }//end for loop j
+                    
+                    //remove from active list 
+                    alignSeqs[i].active = 0;
+                    
+                    out << "ideal_seq_" << (i+1) << '\t' << alignSeqs[i].numIdentical << endl << chunk << endl;;
+                    
+                }//end if active i
+                if(i % 100 == 0)       { m->mothurOutJustToScreen(toString(i) + "\t" + toString(numSeqs - count) + "\t" + toString(count)+"\n");       }
+            }
+        }else {
+            map<int, string> mapFile;
+            map<int, int> originalCount;
+            map<int, int>::iterator itCount;
+            for (int i = 0; i < numSeqs; i++) { mapFile[i] = ""; originalCount[i] = alignSeqs[i].numIdentical; }
+            
+            //think about running through twice...
+            for (int i = 0; i < numSeqs; i++) {
+                
+                //try to merge it into larger seqs
+                for (int j = i+1; j < numSeqs; j++) {
+                    
+                    if (m->control_pressed) { out.close(); return 0; }
+                    
+                    if (originalCount[j] > originalCount[i]) {  //this sequence is more abundant than I am
+                        //are you within "diff" bases
+                        int mismatch = calcMisMatches(alignSeqs[i].seq.getAligned(), alignSeqs[j].seq.getAligned());
+                        
+                        if (mismatch <= diffs) {
+                            //merge
+                            alignSeqs[j].names += ',' + alignSeqs[i].names;
+                            alignSeqs[j].numIdentical += alignSeqs[i].numIdentical;
+                            
+                            mapFile[j] = alignSeqs[i].seq.getName() + "\t" + toString(alignSeqs[i].numIdentical) + "\t" + toString(mismatch) + "\t" + alignSeqs[i].seq.getAligned() + "\n" + mapFile[i];
+                            alignSeqs[i].numIdentical = 0;
+                            originalCount.erase(i);
+                            mapFile[i] = "";
+                            count++;
+                            j+=numSeqs; //exit search, we merged this one in.
+                        }
+                    }//end abundance check
+                }//end for loop j
+                
+                if(i % 100 == 0)       { m->mothurOutJustToScreen(toString(i) + "\t" + toString(numSeqs - count) + "\t" + toString(count)+"\n");       }
+            }
+            
+            for (int i = 0; i < numSeqs; i++) {
+                if (alignSeqs[i].numIdentical != 0) {
+                    out << "ideal_seq_" << (i+1) << '\t' << alignSeqs[i].numIdentical << endl  << alignSeqs[i].seq.getName() + "\t" + toString(alignSeqs[i].numIdentical) + "\t" + toString(0) + "\t" + alignSeqs[i].seq.getAligned() + "\n" << mapFile[i] << endl;
+                }
+            }
+            
+        }
                out.close();
                
                if(numSeqs % 100 != 0)  { m->mothurOut(toString(numSeqs) + "\t" + toString(numSeqs - count) + "\t" + toString(count)); m->mothurOutEndLine();   }       
@@ -861,13 +902,13 @@ void PreClusterCommand::readNameFile(){
                                
                while (!in.eof()) {
                        in >> firstCol >> secondCol; m->gobble(in);
+            
+            m->checkName(firstCol);
+            m->checkName(secondCol);
+            int size = m->getNumNames(secondCol);
+            
                        names[firstCol] = secondCol;
-                       int size = 1;
-
-                       for(int i=0;i<secondCol.size();i++){
-                               if(secondCol[i] == ','){        size++; }
-                       }
-                       sizes[firstCol] = size;
+            sizes[firstCol] = size;
                }
                in.close();
        }