]> git.donarmstrong.com Git - mothur.git/blobdiff - mgclustercommand.cpp
1.12.0
[mothur.git] / mgclustercommand.cpp
index 6d5acb2386ad0ccd836ee41d2d8f5b24c5e4048e..9a1ac215388e09e6477fffa034e0542c38c208b1 100644 (file)
@@ -20,7 +20,7 @@ MGClusterCommand::MGClusterCommand(string option) {
                
                else {
                        //valid paramters for this command
-                       string Array[] =  {"blast", "method", "name", "cutoff", "precision", "length", "min", "penalty", "hcluster","merge","outputdir","inputdir"};
+                       string Array[] =  {"blast", "method", "name", "hard", "cutoff", "precision", "length", "min", "penalty", "hcluster","merge","outputdir","inputdir"};
                        vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
                        
                        OptionParser parser(option);
@@ -104,6 +104,9 @@ MGClusterCommand::MGClusterCommand(string option) {
                        
                        temp = validParameter.validFile(parameters, "hcluster", false);                 if (temp == "not found") { temp = "false"; }
                        hclusterWanted = isTrue(temp); 
+                       
+                       temp = validParameter.validFile(parameters, "hard", false);                     if (temp == "not found") { temp = "F"; }
+                       hard = isTrue(temp);
                }
 
        }
@@ -168,6 +171,8 @@ int MGClusterCommand::execute(){
                
                start = time(NULL);
                oldList = *list;
+               map<string, int> Seq2Bin;
+               map<string, int> oldSeq2Bin;
                
                if (method == "furthest")               { tag = "fn";  }
                else if (method == "nearest")   { tag = "nn";  }
@@ -195,13 +200,15 @@ int MGClusterCommand::execute(){
                        else if(method == "nearest"){   cluster = new SingleLinkage(rabund, list, distMatrix, cutoff, method); }
                        else if(method == "average"){   cluster = new AverageLinkage(rabund, list, distMatrix, cutoff, method); }
                        cluster->setMapWanted(true);
+                       Seq2Bin = cluster->getSeqtoBin();
+                       oldSeq2Bin = Seq2Bin;
                        
                        if (m->control_pressed) { 
                                delete nameMap; delete distMatrix; delete list; delete rabund; delete cluster;
                                listFile.close(); rabundFile.close(); sabundFile.close(); remove((fileroot+ tag + ".list").c_str()); remove((fileroot+ tag + ".rabund").c_str()); remove((fileroot+ tag + ".sabund").c_str());
                                return 0; 
                        }
-                       
+               
                        //cluster using cluster classes
                        while (distMatrix->getSmallDist() < cutoff && distMatrix->getNNodes() > 0){
                                
@@ -214,7 +221,12 @@ int MGClusterCommand::execute(){
                                }
                                
                                float dist = distMatrix->getSmallDist();
-                               float rndDist = roundDist(dist, precision);
+                               float rndDist;
+                               if (hard) {
+                                       rndDist = ceilDist(dist, precision); 
+                               }else{
+                                       rndDist = roundDist(dist, precision); 
+                               }
                                
                                if(previousDist <= 0.0000 && dist != previousDist){
                                        oldList.setLabel("unique");
@@ -222,8 +234,7 @@ int MGClusterCommand::execute(){
                                }
                                else if(rndDist != rndPreviousDist){
                                        if (merge) {
-                                               map<string, int> seq2Bin = cluster->getSeqtoBin();
-                                               ListVector* temp = mergeOPFs(seq2Bin, rndPreviousDist);
+                                               ListVector* temp = mergeOPFs(oldSeq2Bin, rndPreviousDist);
                                                
                                                if (m->control_pressed) { 
                                                        delete nameMap; delete distMatrix; delete list; delete rabund; delete cluster; delete temp;
@@ -239,10 +250,12 @@ int MGClusterCommand::execute(){
                                                printData(&oldList);
                                        }
                                }
-                               
+       
                                previousDist = dist;
                                rndPreviousDist = rndDist;
                                oldList = *list;
+                               Seq2Bin = cluster->getSeqtoBin();
+                               oldSeq2Bin = Seq2Bin;
                        }
                        
                        if(previousDist <= 0.0000){
@@ -251,8 +264,7 @@ int MGClusterCommand::execute(){
                        }
                        else if(rndPreviousDist<cutoff){
                                if (merge) {
-                                       map<string, int> seq2Bin = cluster->getSeqtoBin();
-                                       ListVector* temp = mergeOPFs(seq2Bin, rndPreviousDist);
+                                       ListVector* temp = mergeOPFs(oldSeq2Bin, rndPreviousDist);
                                        
                                        if (m->control_pressed) { 
                                                        delete nameMap; delete distMatrix; delete list; delete rabund; delete cluster; delete temp;
@@ -292,6 +304,8 @@ int MGClusterCommand::execute(){
                        //create cluster
                        hcluster = new HCluster(rabund, list, method, distFile, nameMap, cutoff);
                        hcluster->setMapWanted(true);
+                       Seq2Bin = cluster->getSeqtoBin();
+                       oldSeq2Bin = Seq2Bin;
                        
                        vector<seqDist> seqs; seqs.resize(1); // to start loop
                        //ifstream inHcluster;
@@ -329,7 +343,12 @@ int MGClusterCommand::execute(){
                                                        return 0; 
                                                }
        
-                                               float rndDist = roundDist(seqs[i].dist, precision);
+                                               float rndDist;
+                                               if (hard) {
+                                                       rndDist = ceilDist(seqs[i].dist, precision); 
+                                               }else{
+                                                       rndDist = roundDist(seqs[i].dist, precision); 
+                                               }
                                                                                                
                                                if((previousDist <= 0.0000) && (seqs[i].dist != previousDist)){
                                                        oldList.setLabel("unique");
@@ -337,8 +356,7 @@ int MGClusterCommand::execute(){
                                                }
                                                else if((rndDist != rndPreviousDist)){
                                                        if (merge) {
-                                                               map<string, int> seq2Bin = hcluster->getSeqtoBin();
-                                                               ListVector* temp = mergeOPFs(seq2Bin, rndPreviousDist);
+                                                               ListVector* temp = mergeOPFs(oldSeq2Bin, rndPreviousDist);
                                                                
                                                                if (m->control_pressed) { 
                                                                        delete nameMap;  delete list; delete rabund; delete hcluster; delete temp;
@@ -360,6 +378,8 @@ int MGClusterCommand::execute(){
                                                previousDist = seqs[i].dist;
                                                rndPreviousDist = rndDist;
                                                oldList = *list;
+                                               Seq2Bin = cluster->getSeqtoBin();
+                                               oldSeq2Bin = Seq2Bin;
                                        }
                                }
                        }
@@ -371,8 +391,7 @@ int MGClusterCommand::execute(){
                        }
                        else if(rndPreviousDist<cutoff){
                                if (merge) {
-                                       map<string, int> seq2Bin = hcluster->getSeqtoBin();
-                                       ListVector* temp = mergeOPFs(seq2Bin, rndPreviousDist);
+                                       ListVector* temp = mergeOPFs(oldSeq2Bin, rndPreviousDist);
                                        
                                        if (m->control_pressed) { 
                                                        delete nameMap; delete list; delete rabund; delete hcluster; delete temp;
@@ -452,6 +471,7 @@ ListVector* MGClusterCommand::mergeOPFs(map<string, int> binInfo, float dist){
        try {
                //create new listvector so you don't overwrite the clustering
                ListVector* newList = new ListVector(oldList);
+
                bool done = false;
                ifstream inOverlap;
                int count = 0;
@@ -480,13 +500,16 @@ ListVector* MGClusterCommand::mergeOPFs(map<string, int> binInfo, float dist){
                                        float overlapDistance;
                                        inOverlap >> firstName >> secondName >> overlapDistance; gobble(inOverlap);
                                        
-                                       map<string,int>::iterator itA = nameMap->find(firstName);
-                                       map<string,int>::iterator itB = nameMap->find(secondName);
-                                       if(itA == nameMap->end()){  cerr << "AAError: Sequence '" << firstName << "' was not found in the names file, please correct\n"; exit(1);  }
-                                       if(itB == nameMap->end()){  cerr << "ABError: Sequence '" << secondName << "' was not found in the names file, please correct\n"; exit(1);  }
+                                       //commented out because we check this in readblast already
+                                       //map<string,int>::iterator itA = nameMap->find(firstName);
+                                       //map<string,int>::iterator itB = nameMap->find(secondName);
+                                       //if(itA == nameMap->end()){  cerr << "AAError: Sequence '" << firstName << "' was not found in the names file, please correct\n"; exit(1);  }
+                                       //if(itB == nameMap->end()){  cerr << "ABError: Sequence '" << secondName << "' was not found in the names file, please correct\n"; exit(1);  }
                                        
-                                       overlapNode.seq1 = itA->second;
-                                       overlapNode.seq2 = itB->second;
+                                       //overlapNode.seq1 = itA->second;
+                                       //overlapNode.seq2 = itB->second;
+                                       overlapNode.seq1 = nameMap->get(firstName);
+                                       overlapNode.seq2 = nameMap->get(secondName);
                                        overlapNode.dist = overlapDistance;
                                }else { inOverlap.close(); break; }
                        } 
@@ -495,16 +518,26 @@ ListVector* MGClusterCommand::mergeOPFs(map<string, int> binInfo, float dist){
                                //get names of seqs that overlap
                                string name1 = nameMap->get(overlapNode.seq1);
                                string name2 = nameMap->get(overlapNode.seq2);
-                               
+                       
                                //use binInfo to find out if they are already in the same bin
+                               //map<string, int>::iterator itBin1 = binInfo.find(name1);
+                               //map<string, int>::iterator itBin2 = binInfo.find(name2);
+                               
+                               //if(itBin1 == binInfo.end()){  cerr << "AAError: Sequence '" << name1 << "' does not have any bin info.\n"; exit(1);  }
+                               //if(itBin2 == binInfo.end()){  cerr << "ABError: Sequence '" << name2 << "' does not have any bin info.\n"; exit(1);  }
+
+                               //int binKeep = itBin1->second;
+                               //int binRemove = itBin2->second;
+                               
                                int binKeep = binInfo[name1];
                                int binRemove = binInfo[name2];
-                               
+                       
                                //if not merge bins and update binInfo
                                if(binKeep != binRemove) {
+               
                                        //save names in old bin
-                                       string names = list->get(binRemove);
-                                       
+                                       string names = newList->get(binRemove);
+               
                                        //merge bins into name1s bin
                                        newList->set(binKeep, newList->get(binRemove)+','+newList->get(binKeep));
                                        newList->set(binRemove, "");    
@@ -538,12 +571,12 @@ ListVector* MGClusterCommand::mergeOPFs(map<string, int> binInfo, float dist){
 void MGClusterCommand::sortHclusterFiles(string unsortedDist, string unsortedOverlap) {
        try {
                //sort distFile
-               string sortedDistFile = sortFile(unsortedDist);
+               string sortedDistFile = sortFile(unsortedDist, outputDir);
                remove(unsortedDist.c_str());  //delete unsorted file
                distFile = sortedDistFile;
                
                //sort overlap file
-               string sortedOverlapFile = sortFile(unsortedOverlap);
+               string sortedOverlapFile = sortFile(unsortedOverlap, outputDir);
                remove(unsortedOverlap.c_str());  //delete unsorted file
                overlapFile = sortedOverlapFile;
        }