5 * Created by westcott on 1/7/11.
6 * Copyright 2011 Schloss Lab. All rights reserved.
10 #include "pcacommand.h"
11 #include "inputdata.h"
13 //**********************************************************************************************************************
14 vector<string> PCACommand::getValidParameters(){
16 string Array[] = {"label", "groups","metric","outputdir","inputdir"};
17 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
21 m->errorOut(e, "PCACommand", "getValidParameters");
25 //**********************************************************************************************************************
26 PCACommand::PCACommand(){
28 abort = true; calledHelp = true;
29 vector<string> tempOutNames;
30 outputTypes["pca"] = tempOutNames;
31 outputTypes["loadings"] = tempOutNames;
34 m->errorOut(e, "PCACommand", "PCACommand");
38 //**********************************************************************************************************************
39 vector<string> PCACommand::getRequiredParameters(){
41 vector<string> myArray;
45 m->errorOut(e, "PCACommand", "getRequiredParameters");
49 //**********************************************************************************************************************
50 vector<string> PCACommand::getRequiredFiles(){
52 string Array[] = {"shared","relabund","or"};
53 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
57 m->errorOut(e, "PCACommand", "getRequiredFiles");
61 //**********************************************************************************************************************
63 PCACommand::PCACommand(string option) {
65 abort = false; calledHelp = false;
67 globaldata = GlobalData::getInstance();
69 //allow user to run help
70 if(option == "help") { help(); abort = true; calledHelp = true; }
73 //valid paramters for this command
74 string Array[] = {"label","groups","metric","outputdir", "inputdir"};
75 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
77 OptionParser parser(option);
78 map<string, string> parameters = parser. getParameters();
80 ValidParameters validParameter;
81 map<string, string>::iterator it;
83 //check to make sure all parameters are valid for command
84 for (it = parameters.begin(); it != parameters.end(); it++) {
85 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
87 //if the user changes the input directory command factory will send this info to us in the output parameter
88 string inputDir = validParameter.validFile(parameters, "inputdir", false); if (inputDir == "not found"){ inputDir = ""; }
90 //initialize outputTypes
91 vector<string> tempOutNames;
92 outputTypes["pca"] = tempOutNames;
93 outputTypes["loadings"] = tempOutNames;
95 //make sure the user has already run the read.otu command
96 if ((globaldata->getSharedFile() == "") && (globaldata->getRelAbundFile() == "")) {
97 m->mothurOut("You must read a list and a group, shared or relabund file before you can use the pca command."); m->mothurOutEndLine(); abort = true;
100 if (globaldata->getSharedFile() != "") { mode = "shared"; inputFile = globaldata->getSharedFile(); }
101 if (globaldata->getRelAbundFile() != "") { mode = "relabund"; inputFile = globaldata->getRelAbundFile(); }
103 //if the user changes the output directory command factory will send this info to us in the output parameter
104 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
106 outputDir += m->hasPath(inputFile); //if user entered a file with a path then preserve it
109 string temp = validParameter.validFile(parameters, "metric", false); if (temp == "not found"){ temp = "T"; }
110 metric = m->isTrue(temp);
112 label = validParameter.validFile(parameters, "label", false);
113 if (label == "not found") { label = ""; labels = globaldata->labels; if(labels.size() == 0) { m->mothurOut("You did not provide a label, I will use the first label in your inputfile."); m->mothurOutEndLine(); } }
114 else { m->splitAtDash(label, labels); }
116 groups = validParameter.validFile(parameters, "groups", false);
117 if (groups == "not found") { groups = ""; }
118 else { m->splitAtDash(groups, Groups); }
119 globaldata->Groups = Groups;
124 catch(exception& e) {
125 m->errorOut(e, "PCACommand", "PCACommand");
129 //**********************************************************************************************************************
130 void PCACommand::help(){
132 m->mothurOut("The pca command can only be run after a successful read.otu command of a shared or relabund file."); m->mothurOutEndLine();
133 m->mothurOut("The pca command parameters are label, groups and metric. No parameters are required."); m->mothurOutEndLine();
134 m->mothurOut("The label parameter is used to analyze specific labels in your input. Default is the first label in your shared or relabund file. Multiple labels may be separated by dashes.\n");
135 m->mothurOut("The groups parameter allows you to specify which groups you would like analyzed. Groupnames are separated by dashes.\n");
136 m->mothurOut("The metric parameter allows indicate you if would like the pearson correlation coefficient calculated. Default=True"); m->mothurOutEndLine();
137 m->mothurOut("Example pca(groups=yourGroups).\n");
138 m->mothurOut("Example pca(groups=A-B-C).\n");
139 m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
141 catch(exception& e) {
142 m->errorOut(e, "PCACommand", "help");
146 //**********************************************************************************************************************
147 PCACommand::~PCACommand(){}
148 //**********************************************************************************************************************
149 int PCACommand::execute(){
152 if (abort == true) { if (calledHelp) { return 0; } return 2; }
154 cout.setf(ios::fixed, ios::floatfield);
155 cout.setf(ios::showpoint);
156 cerr.setf(ios::fixed, ios::floatfield);
157 cerr.setf(ios::showpoint);
159 //get first line of shared file
160 vector< vector<double> > matrix;
162 if (mode == "shared") {
163 input = new InputData(inputFile, "sharedfile");
164 }else if (mode == "relabund") {
165 input = new InputData(inputFile, "relabund");
166 }else { m->mothurOut("[ERROR]: filetype not recognized."); m->mothurOutEndLine(); return 0; }
168 vector<SharedRAbundFloatVector*> lookupFloat = input->getSharedRAbundFloatVectors();
169 string lastLabel = lookupFloat[0]->getLabel();
171 set<string> processedLabels;
172 set<string> userLabels = labels;
174 //if the user gave no labels, then use the first one read
175 if (labels.size() == 0) {
178 process(lookupFloat);
181 //as long as you are not at the end of the file or done wih the lines you want
182 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
184 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } delete input; for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } lookupFloat.clear(); return 0; }
186 if(labels.count(lookupFloat[0]->getLabel()) == 1){
187 processedLabels.insert(lookupFloat[0]->getLabel());
188 userLabels.erase(lookupFloat[0]->getLabel());
190 process(lookupFloat);
193 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
194 string saveLabel = lookupFloat[0]->getLabel();
196 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } lookupFloat.clear();
197 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
199 process(lookupFloat);
201 processedLabels.insert(lookupFloat[0]->getLabel());
202 userLabels.erase(lookupFloat[0]->getLabel());
204 //restore real lastlabel to save below
205 lookupFloat[0]->setLabel(saveLabel);
208 lastLabel = lookupFloat[0]->getLabel();
210 //get next line to process
211 //prevent memory leak
212 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } lookupFloat.clear();
213 lookupFloat = input->getSharedRAbundFloatVectors();
217 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } delete input; for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } lookupFloat.clear(); return 0; }
219 //output error messages about any remaining user labels
220 set<string>::iterator it;
221 bool needToRun = false;
222 for (it = userLabels.begin(); it != userLabels.end(); it++) {
223 m->mothurOut("Your file does not include the label " + *it);
224 if (processedLabels.count(lastLabel) != 1) {
225 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
228 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
232 //run last label if you need to
233 if (needToRun == true) {
234 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } } lookupFloat.clear();
235 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
237 process(lookupFloat);
239 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } } lookupFloat.clear();
242 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } } lookupFloat.clear();
245 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
247 m->mothurOutEndLine();
248 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
249 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
250 m->mothurOutEndLine();
254 catch(exception& e) {
255 m->errorOut(e, "PCACommand", "execute");
259 //**********************************************************************************************************************
260 vector< vector<double> > PCACommand::createMatrix(vector<SharedRAbundFloatVector*> lookupFloat){
262 vector< vector<double> > matrix; matrix.resize(lookupFloat.size());
264 //fill matrix with shared files relative abundances
265 for (int i = 0; i < lookupFloat.size(); i++) {
266 for (int j = 0; j < lookupFloat[i]->getNumBins(); j++) {
267 matrix[i].push_back(lookupFloat[i]->getAbundance(j));
271 vector< vector<double> > transposeMatrix; transposeMatrix.resize(matrix[0].size());
272 for (int i = 0; i < transposeMatrix.size(); i++) {
273 for (int j = 0; j < matrix.size(); j++) {
274 transposeMatrix[i].push_back(matrix[j][i]);
278 matrix = linearCalc.matrix_mult(matrix, transposeMatrix);
282 catch(exception& e) {
283 m->errorOut(e, "PCACommand", "createMatrix");
287 //**********************************************************************************************************************
288 int PCACommand::process(vector<SharedRAbundFloatVector*>& lookupFloat){
290 m->mothurOut("\nProcessing " + lookupFloat[0]->getLabel()); m->mothurOutEndLine();
292 vector< vector<double> > matrix; matrix.resize(lookupFloat.size());
294 //fill matrix with shared files relative abundances
295 for (int i = 0; i < lookupFloat.size(); i++) {
296 for (int j = 0; j < lookupFloat[i]->getNumBins(); j++) {
297 matrix[i].push_back(lookupFloat[i]->getAbundance(j));
301 vector< vector<double> > transposeMatrix; transposeMatrix.resize(matrix[0].size());
302 for (int i = 0; i < transposeMatrix.size(); i++) {
303 for (int j = 0; j < matrix.size(); j++) {
304 transposeMatrix[i].push_back(matrix[j][i]);
308 matrix = linearCalc.matrix_mult(matrix, transposeMatrix);
310 double offset = 0.0000;
313 vector<vector<double> > G = matrix;
314 vector<vector<double> > copy_G;
316 for(int count=0;count<2;count++){
317 linearCalc.tred2(G, d, e); if (m->control_pressed) { return 0; }
318 linearCalc.qtli(d, e, G); if (m->control_pressed) { return 0; }
319 offset = d[d.size()-1];
320 if(offset > 0.0) break;
323 if (m->control_pressed) { return 0; }
325 string fbase = outputDir + m->getRootName(m->getSimpleName(inputFile));
326 string outputFileName = fbase + lookupFloat[0]->getLabel();
327 output(outputFileName, globaldata->Groups, G, d);
331 for (int i = 1; i < 4; i++) {
333 vector< vector<double> > EuclidDists = linearCalc.calculateEuclidianDistance(G, i); //G is the pcoa file
335 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
337 double corr = linearCalc.calcPearson(EuclidDists, matrix); //G is the pcoa file, D is the users distance matrix
339 m->mothurOut("Pearson's coefficient using " + toString(i) + " axis: " + toString(corr)); m->mothurOutEndLine();
341 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
347 catch(exception& e) {
348 m->errorOut(e, "PCACommand", "process");
352 /*********************************************************************************************************************************/
354 void PCACommand::output(string fnameRoot, vector<string> name_list, vector<vector<double> >& G, vector<double> d) {
356 int rank = name_list.size();
357 double dsum = 0.0000;
358 for(int i=0;i<rank;i++){
360 for(int j=0;j<rank;j++){
361 if(d[j] >= 0) { G[i][j] *= pow(d[j],0.5); }
362 else { G[i][j] = 0.00000; }
366 ofstream pcaData((fnameRoot+".pca.axes").c_str(), ios::trunc);
367 pcaData.setf(ios::fixed, ios::floatfield);
368 pcaData.setf(ios::showpoint);
369 outputNames.push_back(fnameRoot+".pca.axes");
370 outputTypes["pca"].push_back(fnameRoot+".pca.axes");
372 ofstream pcaLoadings((fnameRoot+".pca.loadings").c_str(), ios::trunc);
373 pcaLoadings.setf(ios::fixed, ios::floatfield);
374 pcaLoadings.setf(ios::showpoint);
375 outputNames.push_back(fnameRoot+".pca.loadings");
376 outputTypes["loadings"].push_back(fnameRoot+".pca.loadings");
378 pcaLoadings << "axis\tloading\n";
379 for(int i=0;i<rank;i++){
380 pcaLoadings << i+1 << '\t' << d[i] * 100.0 / dsum << endl;
384 for(int i=0;i<rank;i++){
385 pcaData << '\t' << "axis" << i+1;
389 for(int i=0;i<rank;i++){
390 pcaData << name_list[i] << '\t';
391 for(int j=0;j<rank;j++){
392 pcaData << G[i][j] << '\t';
397 catch(exception& e) {
398 m->errorOut(e, "PCACommand", "output");
402 /*********************************************************************************************************************************/