2 * matrixoutputcommand.cpp
5 * Created by Sarah Westcott on 5/20/09.
6 * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
10 #include "matrixoutputcommand.h"
11 #include "sharedjabund.h"
12 #include "sharedsorabund.h"
13 #include "sharedjclass.h"
14 #include "sharedsorclass.h"
15 #include "sharedjest.h"
16 #include "sharedsorest.h"
17 #include "sharedthetayc.h"
18 #include "sharedthetan.h"
19 #include "sharedmorisitahorn.h"
20 #include "sharedbraycurtis.h"
23 //**********************************************************************************************************************
25 MatrixOutputCommand::MatrixOutputCommand(string option){
27 globaldata = GlobalData::getInstance();
34 //allow user to run help
35 if(option == "help") { validCalculator = new ValidCalculators(); help(); abort = true; }
38 //valid paramters for this command
39 string Array[] = {"label","calc","groups"};
40 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
42 OptionParser parser(option);
43 map<string,string> parameters = parser.getParameters();
45 ValidParameters validParameter;
47 //check to make sure all parameters are valid for command
48 for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
49 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
52 //make sure the user has already run the read.otu command
53 if (globaldata->getSharedFile() == "") {
54 if (globaldata->getListFile() == "") { mothurOut("You must read a list and a group, or a shared before you can use the dist.shared command."); mothurOutEndLine(); abort = true; }
55 else if (globaldata->getGroupFile() == "") { mothurOut("You must read a list and a group, or a shared before you can use the dist.shared command."); mothurOutEndLine(); abort = true; }
58 //check for optional parameter and set defaults
59 // ...at some point should added some additional type checking...
60 label = validParameter.validFile(parameters, "label", false);
61 if (label == "not found") { label = ""; }
63 if(label != "all") { splitAtDash(label, labels); allLines = 0; }
64 else { allLines = 1; }
67 //if the user has not specified any labels use the ones from read.otu
69 allLines = globaldata->allLines;
70 labels = globaldata->labels;
73 groups = validParameter.validFile(parameters, "groups", false);
74 if (groups == "not found") { groups = ""; }
76 splitAtDash(groups, Groups);
77 globaldata->Groups = Groups;
80 calc = validParameter.validFile(parameters, "calc", false);
81 if (calc == "not found") { calc = "jclass-thetayc"; }
83 if (calc == "default") { calc = "jclass-thetayc"; }
85 splitAtDash(calc, Estimators);
89 validCalculator = new ValidCalculators();
92 for (i=0; i<Estimators.size(); i++) {
93 if (validCalculator->isValidCalculator("matrix", Estimators[i]) == true) {
94 if (Estimators[i] == "jabund") {
95 matrixCalculators.push_back(new JAbund());
96 }else if (Estimators[i] == "sorabund") {
97 matrixCalculators.push_back(new SorAbund());
98 }else if (Estimators[i] == "jclass") {
99 matrixCalculators.push_back(new Jclass());
100 }else if (Estimators[i] == "sorclass") {
101 matrixCalculators.push_back(new SorClass());
102 }else if (Estimators[i] == "jest") {
103 matrixCalculators.push_back(new Jest());
104 }else if (Estimators[i] == "sorest") {
105 matrixCalculators.push_back(new SorEst());
106 }else if (Estimators[i] == "thetayc") {
107 matrixCalculators.push_back(new ThetaYC());
108 }else if (Estimators[i] == "thetan") {
109 matrixCalculators.push_back(new ThetaN());
110 }else if (Estimators[i] == "morisitahorn") {
111 matrixCalculators.push_back(new MorHorn());
112 }else if (Estimators[i] == "braycurtis") {
113 matrixCalculators.push_back(new BrayCurtis());
122 catch(exception& e) {
123 errorOut(e, "MatrixOutputCommand", "MatrixOutputCommand");
128 //**********************************************************************************************************************
130 void MatrixOutputCommand::help(){
132 mothurOut("The dist.shared command can only be executed after a successful read.otu command.\n");
133 mothurOut("The dist.shared command parameters are groups, calc and label. No parameters are required.\n");
134 mothurOut("The groups parameter allows you to specify which of the groups in your groupfile you would like included used.\n");
135 mothurOut("The group names are separated by dashes. The label parameter allows you to select what distance levels you would like distance matrices created for, and is also separated by dashes.\n");
136 mothurOut("The dist.shared command should be in the following format: dist.shared(groups=yourGroups, calc=yourCalcs, label=yourLabels).\n");
137 mothurOut("Example dist.shared(groups=A-B-C, calc=jabund-sorabund).\n");
138 mothurOut("The default value for groups is all the groups in your groupfile.\n");
139 mothurOut("The default value for calc is jclass and thetayc.\n");
140 validCalculator->printCalc("matrix", cout);
141 mothurOut("The dist.shared command outputs a .dist file for each calculator you specify at each distance you choose.\n");
142 mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
144 catch(exception& e) {
145 errorOut(e, "MatrixOutputCommand", "help");
151 //**********************************************************************************************************************
153 MatrixOutputCommand::~MatrixOutputCommand(){
154 if (abort == false) {
155 delete input; globaldata->ginput = NULL;
157 delete validCalculator;
161 //**********************************************************************************************************************
163 int MatrixOutputCommand::execute(){
166 if (abort == true) { return 0; }
168 //if the users entered no valid calculators don't execute command
169 if (matrixCalculators.size() == 0) { mothurOut("No valid calculators."); mothurOutEndLine(); return 0; }
172 read = new ReadOTUFile(globaldata->inputFileName);
173 read->read(&*globaldata);
175 input = globaldata->ginput;
176 lookup = input->getSharedRAbundVectors();
177 string lastLabel = lookup[0]->getLabel();
179 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
180 set<string> processedLabels;
181 set<string> userLabels = labels;
183 if (lookup.size() < 2) { mothurOut("You have not provided enough valid groups. I cannot run the command."); mothurOutEndLine(); return 0;}
185 numGroups = lookup.size();
187 //as long as you are not at the end of the file or done wih the lines you want
188 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
190 if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
191 mothurOut(lookup[0]->getLabel()); mothurOutEndLine();
194 processedLabels.insert(lookup[0]->getLabel());
195 userLabels.erase(lookup[0]->getLabel());
198 if ((anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
199 string saveLabel = lookup[0]->getLabel();
201 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
202 lookup = input->getSharedRAbundVectors(lastLabel);
204 mothurOut(lookup[0]->getLabel()); mothurOutEndLine();
207 processedLabels.insert(lookup[0]->getLabel());
208 userLabels.erase(lookup[0]->getLabel());
210 //restore real lastlabel to save below
211 lookup[0]->setLabel(saveLabel);
214 lastLabel = lookup[0]->getLabel();
216 //get next line to process
217 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
218 lookup = input->getSharedRAbundVectors();
221 //output error messages about any remaining user labels
222 set<string>::iterator it;
223 bool needToRun = false;
224 for (it = userLabels.begin(); it != userLabels.end(); it++) {
225 mothurOut("Your file does not include the label " + *it);
226 if (processedLabels.count(lastLabel) != 1) {
227 mothurOut(". I will use " + lastLabel + "."); mothurOutEndLine();
230 mothurOut(". Please refer to " + lastLabel + "."); mothurOutEndLine();
234 //run last label if you need to
235 if (needToRun == true) {
236 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
237 lookup = input->getSharedRAbundVectors(lastLabel);
239 mothurOut(lookup[0]->getLabel()); mothurOutEndLine();
241 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
245 //reset groups parameter
246 globaldata->Groups.clear();
250 catch(exception& e) {
251 errorOut(e, "MatrixOutputCommand", "execute");
255 /***********************************************************/
256 void MatrixOutputCommand::printSims(ostream& out) {
259 //output column headers
260 out << simMatrix.size() << endl;
262 for (int m = 0; m < simMatrix.size(); m++) {
263 out << lookup[m]->getGroup() << '\t';
264 for (int n = 0; n < m; n++) {
265 out << simMatrix[m][n] << '\t';
271 catch(exception& e) {
272 errorOut(e, "MatrixOutputCommand", "printSims");
276 /***********************************************************/
277 void MatrixOutputCommand::process(vector<SharedRAbundVector*> thisLookup){
281 vector<SharedRAbundVector*> subset;
283 //for each calculator
284 for(int i = 0 ; i < matrixCalculators.size(); i++) {
286 //initialize simMatrix
288 simMatrix.resize(numGroups);
289 for (int m = 0; m < simMatrix.size(); m++) {
290 for (int j = 0; j < simMatrix.size(); j++) {
291 simMatrix[m].push_back(0.0);
295 for (int k = 0; k < thisLookup.size(); k++) {
296 for (int l = k; l < thisLookup.size(); l++) {
297 if (k != l) { //we dont need to similiarity of a groups to itself
298 //get estimated similarity between 2 groups
300 subset.clear(); //clear out old pair of sharedrabunds
301 //add new pair of sharedrabunds
302 subset.push_back(thisLookup[k]); subset.push_back(thisLookup[l]);
304 data = matrixCalculators[i]->getValues(subset); //saves the calculator outputs
305 //save values in similarity matrix
306 simMatrix[k][l] = 1.0 - data[0]; //convert similiarity to distance
307 simMatrix[l][k] = 1.0 - data[0]; //convert similiarity to distance
312 exportFileName = getRootName(globaldata->inputFileName) + matrixCalculators[i]->getName() + "." + thisLookup[0]->getLabel() + ".dist";
313 openOutputFile(exportFileName, out);
322 catch(exception& e) {
323 errorOut(e, "MatrixOutputCommand", "process");
327 /***********************************************************/