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();
35 //allow user to run help
36 if(option == "help") { validCalculator = new ValidCalculators(); help(); abort = true; }
39 //valid paramters for this command
40 string Array[] = {"line","label","calc","groups"};
41 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
43 OptionParser parser(option);
44 map<string,string> parameters = parser.getParameters();
46 ValidParameters validParameter;
48 //check to make sure all parameters are valid for command
49 for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
50 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
53 //make sure the user has already run the read.otu command
54 if (globaldata->getSharedFile() == "") {
55 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; }
56 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; }
59 //check for optional parameter and set defaults
60 // ...at some point should added some additional type checking...
61 line = validParameter.validFile(parameters, "line", false);
62 if (line == "not found") { line = ""; }
64 if(line != "all") { splitAtDash(line, lines); allLines = 0; }
65 else { allLines = 1; }
68 label = validParameter.validFile(parameters, "label", false);
69 if (label == "not found") { label = ""; }
71 if(label != "all") { splitAtDash(label, labels); allLines = 0; }
72 else { allLines = 1; }
75 //make sure user did not use both the line and label parameters
76 if ((line != "") && (label != "")) { mothurOut("You cannot use both the line and label parameters at the same time. "); mothurOutEndLine(); abort = true; }
77 //if the user has not specified any line or labels use the ones from read.otu
78 else if((line == "") && (label == "")) {
79 allLines = globaldata->allLines;
80 labels = globaldata->labels;
81 lines = globaldata->lines;
84 groups = validParameter.validFile(parameters, "groups", false);
85 if (groups == "not found") { groups = ""; }
87 splitAtDash(groups, Groups);
88 globaldata->Groups = Groups;
91 calc = validParameter.validFile(parameters, "calc", false);
92 if (calc == "not found") { calc = "jclass-thetayc"; }
94 if (calc == "default") { calc = "jclass-thetayc"; }
96 splitAtDash(calc, Estimators);
100 validCalculator = new ValidCalculators();
103 for (i=0; i<Estimators.size(); i++) {
104 if (validCalculator->isValidCalculator("matrix", Estimators[i]) == true) {
105 if (Estimators[i] == "jabund") {
106 matrixCalculators.push_back(new JAbund());
107 }else if (Estimators[i] == "sorabund") {
108 matrixCalculators.push_back(new SorAbund());
109 }else if (Estimators[i] == "jclass") {
110 matrixCalculators.push_back(new Jclass());
111 }else if (Estimators[i] == "sorclass") {
112 matrixCalculators.push_back(new SorClass());
113 }else if (Estimators[i] == "jest") {
114 matrixCalculators.push_back(new Jest());
115 }else if (Estimators[i] == "sorest") {
116 matrixCalculators.push_back(new SorEst());
117 }else if (Estimators[i] == "thetayc") {
118 matrixCalculators.push_back(new ThetaYC());
119 }else if (Estimators[i] == "thetan") {
120 matrixCalculators.push_back(new ThetaN());
121 }else if (Estimators[i] == "morisitahorn") {
122 matrixCalculators.push_back(new MorHorn());
123 }else if (Estimators[i] == "braycurtis") {
124 matrixCalculators.push_back(new BrayCurtis());
133 catch(exception& e) {
134 errorOut(e, "MatrixOutputCommand", "MatrixOutputCommand");
139 //**********************************************************************************************************************
141 void MatrixOutputCommand::help(){
143 mothurOut("The dist.shared command can only be executed after a successful read.otu command.\n");
144 mothurOut("The dist.shared command parameters are groups, calc, line and label. The calc parameter is required, and you may not use line and label at the same time.\n");
145 mothurOut("The groups parameter allows you to specify which of the groups in your groupfile you would like included used.\n");
146 mothurOut("The group names are separated by dashes. The line and label allow you to select what distance levels you would like distance matrices created for, and are also separated by dashes.\n");
147 mothurOut("The dist.shared command should be in the following format: dist.shared(groups=yourGroups, calc=yourCalcs, line=yourLines, label=yourLabels).\n");
148 mothurOut("Example dist.shared(groups=A-B-C, line=1-3-5, calc=jabund-sorabund).\n");
149 mothurOut("The default value for groups is all the groups in your groupfile.\n");
150 mothurOut("The default value for calc is jclass and thetayc.\n");
151 validCalculator->printCalc("matrix", cout);
152 mothurOut("The dist.shared command outputs a .dist file for each calculator you specify at each distance you choose.\n");
153 mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
155 catch(exception& e) {
156 errorOut(e, "MatrixOutputCommand", "help");
162 //**********************************************************************************************************************
164 MatrixOutputCommand::~MatrixOutputCommand(){
165 if (abort == false) {
166 delete input; globaldata->ginput = NULL;
168 delete validCalculator;
172 //**********************************************************************************************************************
174 int MatrixOutputCommand::execute(){
177 if (abort == true) { return 0; }
181 //if the users entered no valid calculators don't execute command
182 if (matrixCalculators.size() == 0) { mothurOut("No valid calculators."); mothurOutEndLine(); return 0; }
185 read = new ReadOTUFile(globaldata->inputFileName);
186 read->read(&*globaldata);
188 input = globaldata->ginput;
189 lookup = input->getSharedRAbundVectors();
190 string lastLabel = lookup[0]->getLabel();
192 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
193 set<string> processedLabels;
194 set<string> userLabels = labels;
195 set<int> userLines = lines;
197 if (lookup.size() < 2) { mothurOut("You have not provided enough valid groups. I cannot run the command."); mothurOutEndLine(); return 0;}
199 numGroups = lookup.size();
201 //as long as you are not at the end of the file or done wih the lines you want
202 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0) || (userLines.size() != 0))) {
204 if(allLines == 1 || lines.count(count) == 1 || labels.count(lookup[0]->getLabel()) == 1){
205 mothurOut(lookup[0]->getLabel()); mothurOutEndLine();
208 processedLabels.insert(lookup[0]->getLabel());
209 userLabels.erase(lookup[0]->getLabel());
210 userLines.erase(count);
213 if ((anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
215 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
216 lookup = input->getSharedRAbundVectors(lastLabel);
218 mothurOut(lookup[0]->getLabel()); mothurOutEndLine();
221 processedLabels.insert(lookup[0]->getLabel());
222 userLabels.erase(lookup[0]->getLabel());
226 lastLabel = lookup[0]->getLabel();
228 //get next line to process
229 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
230 lookup = input->getSharedRAbundVectors();
234 //output error messages about any remaining user labels
235 set<string>::iterator it;
236 bool needToRun = false;
237 for (it = userLabels.begin(); it != userLabels.end(); it++) {
238 mothurOut("Your file does not include the label " + *it);
239 if (processedLabels.count(lastLabel) != 1) {
240 mothurOut(". I will use " + lastLabel + "."); mothurOutEndLine();
243 mothurOut(". Please refer to " + lastLabel + "."); mothurOutEndLine();
247 //run last line if you need to
248 if (needToRun == true) {
249 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
250 lookup = input->getSharedRAbundVectors(lastLabel);
252 mothurOut(lookup[0]->getLabel()); mothurOutEndLine();
254 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
258 //reset groups parameter
259 globaldata->Groups.clear();
263 catch(exception& e) {
264 errorOut(e, "MatrixOutputCommand", "execute");
268 /***********************************************************/
269 void MatrixOutputCommand::printSims(ostream& out) {
272 //output column headers
273 out << simMatrix.size() << endl;
275 for (int m = 0; m < simMatrix.size(); m++) {
276 out << lookup[m]->getGroup() << '\t';
277 for (int n = 0; n < m; n++) {
278 out << simMatrix[m][n] << '\t';
284 catch(exception& e) {
285 errorOut(e, "MatrixOutputCommand", "printSims");
289 /***********************************************************/
290 void MatrixOutputCommand::process(vector<SharedRAbundVector*> thisLookup){
294 vector<SharedRAbundVector*> subset;
296 //for each calculator
297 for(int i = 0 ; i < matrixCalculators.size(); i++) {
299 //initialize simMatrix
301 simMatrix.resize(numGroups);
302 for (int m = 0; m < simMatrix.size(); m++) {
303 for (int j = 0; j < simMatrix.size(); j++) {
304 simMatrix[m].push_back(0.0);
308 for (int k = 0; k < thisLookup.size(); k++) {
309 for (int l = k; l < thisLookup.size(); l++) {
310 if (k != l) { //we dont need to similiarity of a groups to itself
311 //get estimated similarity between 2 groups
313 subset.clear(); //clear out old pair of sharedrabunds
314 //add new pair of sharedrabunds
315 subset.push_back(thisLookup[k]); subset.push_back(thisLookup[l]);
317 data = matrixCalculators[i]->getValues(subset); //saves the calculator outputs
318 //save values in similarity matrix
319 simMatrix[k][l] = 1.0 - data[0]; //convert similiarity to distance
320 simMatrix[l][k] = 1.0 - data[0]; //convert similiarity to distance
325 exportFileName = getRootName(globaldata->inputFileName) + matrixCalculators[i]->getName() + "." + thisLookup[0]->getLabel() + ".dist";
326 openOutputFile(exportFileName, out);
335 catch(exception& e) {
336 errorOut(e, "MatrixOutputCommand", "process");
340 /***********************************************************/