]> git.donarmstrong.com Git - mothur.git/blobdiff - unifracunweightedcommand.cpp
merged pat's trim seqs edits with sarah's major overhaul of global data; also added...
[mothur.git] / unifracunweightedcommand.cpp
index 5a7cc367cd12212743f836cda998a10c339b5c09..154c528a851eb297a0783e26bb206aae80fc278c 100644 (file)
 #include "unifracunweightedcommand.h"
 
 /***********************************************************/
-UnifracUnweightedCommand::UnifracUnweightedCommand() {
+UnifracUnweightedCommand::UnifracUnweightedCommand(string option) {
        try {
                globaldata = GlobalData::getInstance();
+               abort = false;
+               Groups.clear();
+               
+               //allow user to run help
+               if(option == "help") { help(); abort = true; }
+               
+               else {
+                       //valid paramters for this command
+                       string Array[] =  {"groups","iters"};
+                       vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+                       
+                       OptionParser parser(option);
+                       map<string,string> parameters = parser.getParameters();
+                       
+                       ValidParameters validParameter;
+               
+                       //check to make sure all parameters are valid for command
+                       for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) { 
+                               if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
+                       }
+                       
+                       if (globaldata->gTree.size() == 0) {//no trees were read
+                               cout << "You must execute the read.tree command, before you may execute the unifrac.unweighted command." << endl; abort = true;  }
+                                                                               
+                       //check for optional parameter and set defaults
+                       // ...at some point should added some additional type checking...
+                       groups = validParameter.validFile(parameters, "groups", false);                 
+                       if (groups == "not found") { groups = ""; }
+                       else { 
+                               splitAtDash(groups, Groups);
+                               globaldata->Groups = Groups;
+                       }
+                               
+                       itersString = validParameter.validFile(parameters, "iters", false);                     if (itersString == "not found") { itersString = "1000"; }
+                       convert(itersString, iters); 
+                       
+                       
+                       if (abort == false) {
+                               T = globaldata->gTree;
+                               tmap = globaldata->gTreemap;
+                               sumFile = globaldata->getTreeFile() + ".uwsummary";
+                               openOutputFile(sumFile, outSum);
+                               
+                               util = new SharedUtil();
+                               util->setGroups(globaldata->Groups, tmap->namesOfGroups, allGroups, numGroups, "unweighted");   //sets the groups the user wants to analyze
+                               util->getCombos(groupComb, globaldata->Groups, numComp);
+                               
+                               if (numGroups == 1) { numComp++; groupComb.push_back(allGroups); }
+                               
+                               unweighted = new Unweighted(tmap);
+                               
+                       }
+                       
+               }
                
-               T = globaldata->gTree;
-               tmap = globaldata->gTreemap;
-               sumFile = globaldata->getTreeFile() + ".uwsummary";
-               openOutputFile(sumFile, outSum);
-
-               setGroups(); //sets users groups to analyze
-               convert(globaldata->getIters(), iters);  //how many random trees to generate
-               unweighted = new Unweighted(tmap);
-
        }
        catch(exception& e) {
                cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function UnifracUnweightedCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
@@ -33,39 +78,64 @@ UnifracUnweightedCommand::UnifracUnweightedCommand() {
                exit(1);
        }
 }
+
+//**********************************************************************************************************************
+
+void UnifracUnweightedCommand::help(){
+       try {
+               cout << "The unifrac.unweighted command can only be executed after a successful read.tree command." << "\n";
+               cout << "The unifrac.unweighted command parameters are groups and iters.  No parameters are required." << "\n";
+               cout << "The groups parameter allows you to specify which of the groups in your groupfile you would like analyzed.  You must enter at least 1 valid group." << "\n";
+               cout << "The group names are separated by dashes.  The iters parameter allows you to specify how many random trees you would like compared to your tree." << "\n";
+               cout << "The unifrac.unweighted command should be in the following format: unifrac.unweighted(groups=yourGroups, iters=yourIters)." << "\n";
+               cout << "Example unifrac.unweighted(groups=A-B-C, iters=500)." << "\n";
+               cout << "The default value for groups is all the groups in your groupfile, and iters is 1000." << "\n";
+               cout << "The unifrac.unweighted command output two files: .unweighted and .uwsummary their descriptions are in the manual." << "\n";
+               cout << "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups)." << "\n" << "\n";
+       }
+       catch(exception& e) {
+               cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function help. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }
+       catch(...) {
+               cout << "An unknown error has occurred in the UnifracUnweightedCommand class function help. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               exit(1);
+       }       
+}
+
+
 /***********************************************************/
 int UnifracUnweightedCommand::execute() {
        try {
-       
+               
+               if (abort == true) { return 0; }
+               
                userData.resize(numComp,0);  //data[0] = unweightedscore 
                randomData.resize(numComp,0); //data[0] = unweightedscore
                //create new tree with same num nodes and leaves as users
-                               
+               
+               outSum << "Tree#" << '\t' << "Groups" << '\t'  <<  "UWScore" <<'\t' << "UWSig" <<  endl;
+               cout << "Tree#" << '\t' << "Groups" << '\t'  <<  "UWScore" << '\t' << "UWSig" <<  endl;
+               
                //get pscores for users trees
                for (int i = 0; i < T.size(); i++) {
-                       unweightedFile = globaldata->getTreeFile()  + toString(i+1) + ".unweighted";
-                       unweightedFileout = globaldata->getTreeFile() + "temp." + toString(i+1) + ".unweighted";
-                                               //column headers
-                       outSum << "Tree# " << i+1 << endl;
-                       outSum << "Comb" << '\t'  <<  "UWScore" << '\t' << '\t' << "UWSig" <<  endl;
-
+                       counter = 0;
+                       
+                       output = new ColumnFile(globaldata->getTreeFile()  + toString(i+1) + ".unweighted", itersString);
+                       
                        //get unweighted for users tree
                        rscoreFreq.resize(numComp);  
                        rCumul.resize(numComp);  
-                       validScores.resize(numComp); 
                        utreeScores.resize(numComp);  
                        UWScoreSig.resize(numComp); 
 
-                       cout << "Processing tree " << i+1 << endl;
                        userData = unweighted->getValues(T[i]);  //userData[0] = unweightedscore
                        
                        //output scores for each combination
                        for(int k = 0; k < numComp; k++) {
-                               //add users score to valid scores
-                               validScores[k][userData[k]] = userData[k];
-                       
                                //saves users score
                                utreeScores[k].push_back(userData[k]);
+
                        }
                        
                        //get unweighted scores for random trees
@@ -74,45 +144,49 @@ int UnifracUnweightedCommand::execute() {
                                randomData = unweighted->getValues(T[i], "", "");
                                
                                for(int k = 0; k < numComp; k++) {      
-//cout << "iter " << j << " comp " << k << " = " << randomData[k] << endl;
                                        //add trees unweighted score to map of scores
-                                       it2 = rscoreFreq[k].find(randomData[k]);
-                                       if (it2 != rscoreFreq[k].end()) {//already have that score
+                                       map<float,float>::iterator it = rscoreFreq[k].find(randomData[k]);
+                                       if (it != rscoreFreq[k].end()) {//already have that score
                                                rscoreFreq[k][randomData[k]]++;
                                        }else{//first time we have seen this score
                                                rscoreFreq[k][randomData[k]] = 1;
                                        }
                                
                                        //add randoms score to validscores
-                                       validScores[k][randomData[k]] = randomData[k];
+                                       validScores[randomData[k]] = randomData[k];
                                }
+                               
                        }
                
-               for(int a = 0; a < numComp; a++) {
-                       float rcumul = 1.0000;
-                       //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
-                       for (it = validScores[a].begin(); it != validScores[a].end(); it++) { 
-                               //make rscoreFreq map and rCumul
-                               it2 = rscoreFreq[a].find(it->first);
-                               rCumul[a][it->first] = rcumul;
-                               //get percentage of random trees with that info
-                               if (it2 != rscoreFreq[a].end()) {  rscoreFreq[a][it->first] /= iters; rcumul-= it2->second;  }
-                               else { rscoreFreq[a][it->first] = 0.0000; } //no random trees with that score
+                       for(int a = 0; a < numComp; a++) {
+                               float rcumul = 1.0000;
+                               //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
+                               for (map<float,float>::iterator it = validScores.begin(); it != validScores.end(); it++) { 
+                                       //make rscoreFreq map and rCumul
+                                       map<float,float>::iterator it2 = rscoreFreq[a].find(it->first);
+                                       rCumul[a][it->first] = rcumul;
+                                       //get percentage of random trees with that info
+                                       if (it2 != rscoreFreq[a].end()) {  rscoreFreq[a][it->first] /= iters; rcumul-= it2->second;  }
+                                       else { rscoreFreq[a][it->first] = 0.0000; } //no random trees with that score
+                               }
+                               UWScoreSig[a].push_back(rCumul[a][userData[a]]);
                        }
-                       UWScoreSig[a].push_back(rCumul[a][userData[a]]);
-               }
                
-               printUnweightedFile();
-               printUWSummaryFile();
                
-               rscoreFreq.clear(); 
-               rCumul.clear();  
-               validScores.clear(); 
-               utreeScores.clear();  
-               UWScoreSig.clear(); 
-       }
+               
+                       printUnweightedFile();
+                       printUWSummaryFile(i);
+                       
+                       delete output;
+                       rscoreFreq.clear(); 
+                       rCumul.clear();  
+                       validScores.clear(); 
+                       utreeScores.clear();  
+                       UWScoreSig.clear(); 
+               }
+               
                //reset groups parameter
-               globaldata->Groups.clear(); globaldata->setGroups("");
+               globaldata->Groups.clear(); 
                outSum.close();
                
                return 0;
@@ -131,22 +205,19 @@ int UnifracUnweightedCommand::execute() {
 void UnifracUnweightedCommand::printUnweightedFile() {
        try {
                vector<double> data;
+               vector<string> tags;
+               tags.push_back("Score"); tags.push_back("RandFreq"); tags.push_back("RandCumul");
                
                for(int a = 0; a < numComp; a++) {
-                       initFile(groupComb[a]);
+                       output->initFile(groupComb[a], tags);
                        //print each line
-                       for (it = validScores[a].begin(); it != validScores[a].end(); it++) { 
+                       for (map<float,float>::iterator it = validScores.begin(); it != validScores.end(); it++) { 
                                data.push_back(it->first);  data.push_back(rscoreFreq[a][it->first]); data.push_back(rCumul[a][it->first]); 
-                               output(data);
+                               output->output(data);
                                data.clear();
                        } 
-                       resetFile();
+                       output->resetFile();
                }
-               
-               out.close();
-               inFile.close();
-               remove(unweightedFileout.c_str());
-               
        }
        catch(exception& e) {
                cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function printUnweightedFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
@@ -159,178 +230,38 @@ void UnifracUnweightedCommand::printUnweightedFile() {
 }
 
 /***********************************************************/
-void UnifracUnweightedCommand::printUWSummaryFile() {
+void UnifracUnweightedCommand::printUWSummaryFile(int i) {
        try {
                                
                //format output
                outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
-               
+                       
                //print each line
-               for(int a = 0; a < numComp; a++) {
-                       outSum << setprecision(6) << groupComb[a] << '\t' << '\t' << utreeScores[a][0] << '\t' << UWScoreSig[a][0] << endl;
-                       cout << setprecision(6)  << groupComb[a] << '\t' << '\t' << utreeScores[a][0] << '\t' << UWScoreSig[a][0] << endl; 
-               }
-               
-       }
-       catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function printUWSummaryFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }
-       catch(...) {
-               cout << "An unknown error has occurred in the UnifracUnweightedCommand class function printUWSummaryFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }
-}
-/***********************************************************/
 
-void UnifracUnweightedCommand::setGroups() {
-       try {
-               string allGroups = "";
-               numGroups = 0;
-               //if the user has not entered specific groups to analyze then do them all
-               if (globaldata->Groups.size() != 0) {
-                       if (globaldata->Groups[0] != "all") {
-                               //check that groups are valid
-                               for (int i = 0; i < globaldata->Groups.size(); i++) {
-                                       if (tmap->isValidGroup(globaldata->Groups[i]) != true) {
-                                               cout << globaldata->Groups[i] << " is not a valid group, and will be disregarded." << endl;
-                                               // erase the invalid group from globaldata->Groups
-                                               globaldata->Groups.erase(globaldata->Groups.begin()+i);
-                                       }
-                               }
+               for(int a = 0; a < numComp; a++) {
+                       outSum << i+1 << '\t';
+                       cout << i+1 << '\t';
                        
-                               //if the user only entered invalid groups
-                               if (globaldata->Groups.size() == 0) { 
-                                       cout << "When using the groups parameter you must have at least 1 valid group. I will run the command using all the groups in your groupfile." << endl; 
-                                       for (int i = 0; i < tmap->namesOfGroups.size(); i++) {
-                                               globaldata->Groups.push_back(tmap->namesOfGroups[i]);
-                                               numGroups++;
-                                               allGroups += tmap->namesOfGroups[i];
-                                       }
-                               }else {
-                                       for (int i = 0; i < globaldata->Groups.size(); i++) {
-                                               allGroups += globaldata->Groups[i];
-                                               numGroups++;
-                                       }
-                               }
-                       }else{//user has enter "all" and wants the default groups
-                               for (int i = 0; i < tmap->namesOfGroups.size(); i++) {
-                                       globaldata->Groups.push_back(tmap->namesOfGroups[i]);
-                                       numGroups++;
-                                       allGroups += tmap->namesOfGroups[i];
-                               }
-                               globaldata->setGroups("");
-                       }
-               }else {
-                       for (int i = 0; i < tmap->namesOfGroups.size(); i++) {
-                               allGroups += tmap->namesOfGroups[i];
+                       if (UWScoreSig[a][0] > (1/(float)iters)) {
+                               outSum << setprecision(6) << groupComb[a]  << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << UWScoreSig[a][0] << endl;
+                               cout << setprecision(6)  << groupComb[a]  << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << UWScoreSig[a][0] << endl; 
+                       }else {
+                               outSum << setprecision(6) << groupComb[a]  << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << "<" << (1/float(iters)) << endl;
+                               cout << setprecision(6)  << groupComb[a]  << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << "<" << (1/float(iters)) << endl; 
                        }
-                       numGroups = 1;
                }
                
-               //calculate number of comparsions
-               numComp = 0;
-               for (int r=0; r<numGroups; r++) { 
-                       for (int l = r+1; l < numGroups; l++) {
-                               groupComb.push_back(globaldata->Groups[r]+globaldata->Groups[l]);
-                               numComp++;
-                       }
-               }
-               
-               //ABC
-               if (numComp != 1) {
-                       groupComb.push_back(allGroups);
-                       numComp++;
-               }
        }
        catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function setGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }
-       catch(...) {
-               cout << "An unknown error has occurred in the UnifracUnweightedCommand class function setGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }               
-
-}
-/*****************************************************************/
-
-void UnifracUnweightedCommand::initFile(string label){
-       try {
-               if(counter != 0){
-                       openOutputFile(unweightedFileout, out);
-                       openInputFile(unweightedFile, inFile);
-
-                       string inputBuffer;
-                       getline(inFile, inputBuffer);
-               
-                       out     <<  inputBuffer << '\t' << label + "Score" << '\t' << label + "RandFreq" << '\t' << label + "RandCumul" << endl;                
-               }else{
-                       openOutputFile(unweightedFileout, out);
-                       out     << label + "Score" << '\t' << label + "RandFreq" << '\t' << label + "RandCumul" << endl;
-               }
-
-               out.setf(ios::fixed, ios::floatfield);
-               out.setf(ios::showpoint);
-       }
-       catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function initFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function printUWSummaryFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
                exit(1);
        }
        catch(...) {
-               cout << "An unknown error has occurred in the UnifracUnweightedCommand class function initFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               cout << "An unknown error has occurred in the UnifracUnweightedCommand class function printUWSummaryFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
                exit(1);
        }
 }
 
-/***********************************************************************/
-
-void UnifracUnweightedCommand::output(vector<double> data){
-       try {
-               if(counter != 0){               
-                       string inputBuffer;
-                       getline(inFile, inputBuffer);
-               
-                       out     <<  inputBuffer << setprecision(6) << '\t' << data[0] << '\t' << data[1] << '\t' << data[2] << endl;
-               }
-               else{
-                       out << setprecision(6) << data[0] << '\t' << data[1] << '\t' << data[2] << endl;
-               }
-
-       }
-       catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function output. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }
-       catch(...) {
-               cout << "An unknown error has occurred in the UnifracUnweightedCommand class function output. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }
-};
-
-/***********************************************************************/
+/***********************************************************/
 
-void UnifracUnweightedCommand::resetFile(){
-       try {
-               if(counter != 0){
-                       out.close();
-                       inFile.close();
-               }
-               else{
-                       out.close();
-               }
-               counter = 1;
-               
-               remove(unweightedFile.c_str());
-               rename(unweightedFileout.c_str(), unweightedFile.c_str());
-       }
-       catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the UnifracUnweightedCommand class Function resetFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }
-       catch(...) {
-               cout << "An unknown error has occurred in the UnifracUnweightedCommand class function resetFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }       
-}