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","outputdir","inputdir", "output"};
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 //if the user changes the output directory command factory will send this info to us in the output parameter
53 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
55 outputDir += hasPath(globaldata->inputFileName); //if user entered a file with a path then preserve it
58 //make sure the user has already run the read.otu command
59 if (globaldata->getSharedFile() == "") {
60 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; }
61 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; }
64 //check for optional parameter and set defaults
65 // ...at some point should added some additional type checking...
66 label = validParameter.validFile(parameters, "label", false);
67 if (label == "not found") { label = ""; }
69 if(label != "all") { splitAtDash(label, labels); allLines = 0; }
70 else { allLines = 1; }
73 output = validParameter.validFile(parameters, "output", false); if(output == "not found"){ output = "lt"; }
74 if ((output != "lt") && (output != "square")) { mothurOut(output + " is not a valid output form. Options are lt and square. I will use lt."); mothurOutEndLine(); output = "lt"; }
76 //if the user has not specified any labels use the ones from read.otu
78 allLines = globaldata->allLines;
79 labels = globaldata->labels;
82 groups = validParameter.validFile(parameters, "groups", false);
83 if (groups == "not found") { groups = ""; }
85 splitAtDash(groups, Groups);
86 globaldata->Groups = Groups;
89 calc = validParameter.validFile(parameters, "calc", false);
90 if (calc == "not found") { calc = "jclass-thetayc"; }
92 if (calc == "default") { calc = "jclass-thetayc"; }
94 splitAtDash(calc, Estimators);
98 validCalculator = new ValidCalculators();
101 for (i=0; i<Estimators.size(); i++) {
102 if (validCalculator->isValidCalculator("matrix", Estimators[i]) == true) {
103 if (Estimators[i] == "jabund") {
104 matrixCalculators.push_back(new JAbund());
105 }else if (Estimators[i] == "sorabund") {
106 matrixCalculators.push_back(new SorAbund());
107 }else if (Estimators[i] == "jclass") {
108 matrixCalculators.push_back(new Jclass());
109 }else if (Estimators[i] == "sorclass") {
110 matrixCalculators.push_back(new SorClass());
111 }else if (Estimators[i] == "jest") {
112 matrixCalculators.push_back(new Jest());
113 }else if (Estimators[i] == "sorest") {
114 matrixCalculators.push_back(new SorEst());
115 }else if (Estimators[i] == "thetayc") {
116 matrixCalculators.push_back(new ThetaYC());
117 }else if (Estimators[i] == "thetan") {
118 matrixCalculators.push_back(new ThetaN());
119 }else if (Estimators[i] == "morisitahorn") {
120 matrixCalculators.push_back(new MorHorn());
121 }else if (Estimators[i] == "braycurtis") {
122 matrixCalculators.push_back(new BrayCurtis());
131 catch(exception& e) {
132 errorOut(e, "MatrixOutputCommand", "MatrixOutputCommand");
137 //**********************************************************************************************************************
139 void MatrixOutputCommand::help(){
141 mothurOut("The dist.shared command can only be executed after a successful read.otu command.\n");
142 mothurOut("The dist.shared command parameters are groups, calc, output and label. No parameters are required.\n");
143 mothurOut("The groups parameter allows you to specify which of the groups in your groupfile you would like included used.\n");
144 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");
145 mothurOut("The dist.shared command should be in the following format: dist.shared(groups=yourGroups, calc=yourCalcs, label=yourLabels).\n");
146 mothurOut("The output parameter allows you to specify format of your distance matrix. Options are lt, and square. The default is lt.\n");
147 mothurOut("Example dist.shared(groups=A-B-C, calc=jabund-sorabund).\n");
148 mothurOut("The default value for groups is all the groups in your groupfile.\n");
149 mothurOut("The default value for calc is jclass and thetayc.\n");
150 validCalculator->printCalc("matrix", cout);
151 mothurOut("The dist.shared command outputs a .dist file for each calculator you specify at each distance you choose.\n");
152 mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
154 catch(exception& e) {
155 errorOut(e, "MatrixOutputCommand", "help");
161 //**********************************************************************************************************************
163 MatrixOutputCommand::~MatrixOutputCommand(){
164 if (abort == false) {
165 delete input; globaldata->ginput = NULL;
167 delete validCalculator;
171 //**********************************************************************************************************************
173 int MatrixOutputCommand::execute(){
176 if (abort == true) { return 0; }
178 //if the users entered no valid calculators don't execute command
179 if (matrixCalculators.size() == 0) { mothurOut("No valid calculators."); mothurOutEndLine(); return 0; }
182 read = new ReadOTUFile(globaldata->inputFileName);
183 read->read(&*globaldata);
185 input = globaldata->ginput;
186 lookup = input->getSharedRAbundVectors();
187 string lastLabel = lookup[0]->getLabel();
189 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
190 set<string> processedLabels;
191 set<string> userLabels = labels;
193 if (lookup.size() < 2) { mothurOut("You have not provided enough valid groups. I cannot run the command."); mothurOutEndLine(); return 0;}
195 numGroups = lookup.size();
197 //as long as you are not at the end of the file or done wih the lines you want
198 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
200 if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
201 mothurOut(lookup[0]->getLabel()); mothurOutEndLine();
204 processedLabels.insert(lookup[0]->getLabel());
205 userLabels.erase(lookup[0]->getLabel());
208 if ((anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
209 string saveLabel = lookup[0]->getLabel();
211 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
212 lookup = input->getSharedRAbundVectors(lastLabel);
214 mothurOut(lookup[0]->getLabel()); mothurOutEndLine();
217 processedLabels.insert(lookup[0]->getLabel());
218 userLabels.erase(lookup[0]->getLabel());
220 //restore real lastlabel to save below
221 lookup[0]->setLabel(saveLabel);
224 lastLabel = lookup[0]->getLabel();
226 //get next line to process
227 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
228 lookup = input->getSharedRAbundVectors();
231 //output error messages about any remaining user labels
232 set<string>::iterator it;
233 bool needToRun = false;
234 for (it = userLabels.begin(); it != userLabels.end(); it++) {
235 mothurOut("Your file does not include the label " + *it);
236 if (processedLabels.count(lastLabel) != 1) {
237 mothurOut(". I will use " + lastLabel + "."); mothurOutEndLine();
240 mothurOut(". Please refer to " + lastLabel + "."); mothurOutEndLine();
244 //run last label if you need to
245 if (needToRun == true) {
246 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
247 lookup = input->getSharedRAbundVectors(lastLabel);
249 mothurOut(lookup[0]->getLabel()); mothurOutEndLine();
251 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
255 //reset groups parameter
256 globaldata->Groups.clear();
260 catch(exception& e) {
261 errorOut(e, "MatrixOutputCommand", "execute");
265 /***********************************************************/
266 void MatrixOutputCommand::printSims(ostream& out) {
269 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
272 out << simMatrix.size() << endl;
274 if (output == "lt") {
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';
283 for (int m = 0; m < simMatrix.size(); m++) {
284 out << lookup[m]->getGroup() << '\t';
285 for (int n = 0; n < simMatrix[m].size(); n++) {
286 out << simMatrix[m][n] << '\t';
292 catch(exception& e) {
293 errorOut(e, "MatrixOutputCommand", "printSims");
297 /***********************************************************/
298 void MatrixOutputCommand::process(vector<SharedRAbundVector*> thisLookup){
302 vector<SharedRAbundVector*> subset;
304 //for each calculator
305 for(int i = 0 ; i < matrixCalculators.size(); i++) {
307 //initialize simMatrix
309 simMatrix.resize(numGroups);
310 for (int m = 0; m < simMatrix.size(); m++) {
311 for (int j = 0; j < simMatrix.size(); j++) {
312 simMatrix[m].push_back(0.0);
316 for (int k = 0; k < thisLookup.size(); k++) {
317 for (int l = k; l < thisLookup.size(); l++) {
318 if (k != l) { //we dont need to similiarity of a groups to itself
319 //get estimated similarity between 2 groups
321 subset.clear(); //clear out old pair of sharedrabunds
322 //add new pair of sharedrabunds
323 subset.push_back(thisLookup[k]); subset.push_back(thisLookup[l]);
325 data = matrixCalculators[i]->getValues(subset); //saves the calculator outputs
326 //save values in similarity matrix
327 simMatrix[k][l] = 1.0 - data[0]; //convert similiarity to distance
328 simMatrix[l][k] = 1.0 - data[0]; //convert similiarity to distance
333 exportFileName = outputDir + getRootName(getSimpleName(globaldata->inputFileName)) + matrixCalculators[i]->getName() + "." + thisLookup[0]->getLabel() + "." + output + ".dist";
334 openOutputFile(exportFileName, out);
343 catch(exception& e) {
344 errorOut(e, "MatrixOutputCommand", "process");
348 /***********************************************************/