1
2 /*
3 * pcacommand.cpp
4 * Mothur
5 *
6 * Created by westcott on 1/4/10.
7 * Copyright 2010 Schloss Lab. All rights reserved.
8 *
9 */
10
11 #include "pcoacommand.h"
12 #include "readphylipvector.h"
13
14
15 //**********************************************************************************************************************
setParameters()16 vector<string> PCOACommand::setParameters(){
17 try {
18 CommandParameter pphylip("phylip", "InputTypes", "", "", "none", "none", "none","pcoa-loadings",false,true,true); parameters.push_back(pphylip);
19 CommandParameter pmetric("metric", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(pmetric);
20 CommandParameter pseed("seed", "Number", "", "0", "", "", "","",false,false); parameters.push_back(pseed);
21 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
22 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
23
24 abort = false; calledHelp = false;
25
26 vector<string> tempOutNames;
27 outputTypes["pcoa"] = tempOutNames;
28 outputTypes["loadings"] = tempOutNames;
29
30 vector<string> myArray;
31 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
32 return myArray;
33 }
34 catch(exception& e) {
35 m->errorOut(e, "PCOACommand", "setParameters");
36 exit(1);
37 }
38 }
39 //**********************************************************************************************************************
getHelpString()40 string PCOACommand::getHelpString(){
41 try {
42 string helpString = "";
43 helpString += "The pcoa command parameters are phylip and metric";
44 helpString += "The phylip parameter allows you to enter your distance file.";
45 helpString += "The metric parameter allows indicate you if would like the pearson correlation coefficient calculated. Default=True";
46 helpString += "Example pcoa(phylip=yourDistanceFile).\n";
47 return helpString;
48 }
49 catch(exception& e) {
50 m->errorOut(e, "PCOACommand", "getHelpString");
51 exit(1);
52 }
53 }
54 //**********************************************************************************************************************
getOutputPattern(string type)55 string PCOACommand::getOutputPattern(string type) {
56 try {
57 string pattern = "";
58
59 if (type == "pcoa") { pattern = "[filename],pcoa.axes"; }
60 else if (type == "loadings") { pattern = "[filename],pcoa.loadings"; }
61 else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->setControl_pressed(true); }
62
63 return pattern;
64 }
65 catch(exception& e) {
66 m->errorOut(e, "PCOACommand", "getOutputPattern");
67 exit(1);
68 }
69 }
70
71 //**********************************************************************************************************************
PCOACommand(string option)72 PCOACommand::PCOACommand(string option) : Command() {
73 try {
74 if(option == "help") { help(); abort = true; calledHelp = true; }
75 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
76 else if(option == "category") { abort = true; calledHelp = true; }
77
78 else {
79 OptionParser parser(option, setParameters());
80 map<string, string> parameters = parser. getParameters();
81
82 ValidParameters validParameter;
83 phylipfile = validParameter.validFile(parameters, "phylip");
84 if (phylipfile == "not open") { abort = true; }
85 else if (phylipfile == "not found") {
86 //if there is a current phylip file, use it
87 phylipfile = current->getPhylipFile();
88 if (phylipfile != "") { m->mothurOut("Using " + phylipfile + " as input file for the phylip parameter.\n"); }
89 else { m->mothurOut("You have no current phylip file and the phylip parameter is required.\n"); abort = true; }
90 }else { current->setPhylipFile(phylipfile); }
91
92 filename = phylipfile;
93
94 if (outputdir == ""){ outputdir += util.hasPath(phylipfile); }
95
96 string temp = validParameter.valid(parameters, "metric"); if (temp == "not found"){ temp = "T"; }
97 metric = util.isTrue(temp);
98 }
99
100 }
101 catch(exception& e) {
102 m->errorOut(e, "PCOACommand", "PCOACommand");
103 exit(1);
104 }
105 }
106 //**********************************************************************************************************************
execute()107 int PCOACommand::execute(){
108 try {
109
110 if (abort) { if (calledHelp) { return 0; } return 2; }
111
112 cout.setf(ios::fixed, ios::floatfield);
113 cout.setf(ios::showpoint);
114 cerr.setf(ios::fixed, ios::floatfield);
115 cerr.setf(ios::showpoint);
116
117 vector<string> names;
118 vector<vector<double> > D;
119
120 fbase = outputdir + util.getRootName(util.getSimpleName(filename));
121
122 ReadPhylipVector readFile(filename);
123 names = readFile.read(D);
124
125 if (m->getControl_pressed()) { return 0; }
126
127 double offset = 0.0000;
128 vector<double> d;
129 vector<double> e;
130 vector<vector<double> > G = D;
131 //vector<vector<double> > copy_G;
132
133 m->mothurOut("\nProcessing...\n");
134
135 for(int count=0;count<2;count++){
136 linearCalc.recenter(offset, D, G); if (m->getControl_pressed()) { return 0; }
137 linearCalc.tred2(G, d, e); if (m->getControl_pressed()) { return 0; }
138 linearCalc.qtli(d, e, G); if (m->getControl_pressed()) { return 0; }
139 offset = d[d.size()-1];
140 if(offset > 0.0) break;
141 }
142
143 if (m->getControl_pressed()) { return 0; }
144
145 output(fbase, names, G, d);
146
147 if (m->getControl_pressed()) { for (int i = 0; i < outputNames.size(); i++) { util.mothurRemove(outputNames[i]); } return 0; }
148
149 if (metric) {
150
151 for (int i = 1; i < 4; i++) {
152
153 vector< vector<double> > EuclidDists = linearCalc.calculateEuclidianDistance(G, i); //G is the pcoa file
154
155 if (m->getControl_pressed()) { for (int i = 0; i < outputNames.size(); i++) { util.mothurRemove(outputNames[i]); } return 0; }
156
157 double corr = linearCalc.calcPearson(EuclidDists, D); //G is the pcoa file, D is the users distance matrix
158
159 m->mothurOut("Rsq " + toString(i) + " axis: " + toString(corr * corr)); m->mothurOutEndLine();
160
161 if (m->getControl_pressed()) { for (int i = 0; i < outputNames.size(); i++) { util.mothurRemove(outputNames[i]); } return 0; }
162 }
163 }
164
165 m->mothurOut("\nOutput File Names: \n");
166 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i] +"\n"); } m->mothurOutEndLine();
167
168 return 0;
169 }
170 catch(exception& e) {
171 m->errorOut(e, "PCOACommand", "execute");
172 exit(1);
173 }
174 }
175 /*********************************************************************************************************************************/
176
get_comment(istream & f,char begin,char end)177 void PCOACommand::get_comment(istream& f, char begin, char end){
178 try {
179 char d=f.get();
180 while(d != end){ d = f.get(); }
181 d = f.peek();
182 }
183 catch(exception& e) {
184 m->errorOut(e, "PCOACommand", "get_comment");
185 exit(1);
186 }
187 }
188 /*********************************************************************************************************************************/
189
output(string fnameRoot,vector<string> name_list,vector<vector<double>> & G,vector<double> d)190 void PCOACommand::output(string fnameRoot, vector<string> name_list, vector<vector<double> >& G, vector<double> d) {
191 try {
192 int rank = name_list.size();
193 double dsum = 0.0000;
194 for(int i=0;i<rank;i++){
195 dsum += d[i];
196 for(int j=0;j<rank;j++){
197 if(d[j] >= 0) { G[i][j] *= pow(d[j],0.5); }
198 else { G[i][j] = 0.00000; }
199 }
200 }
201
202 ofstream pcaData;
203 map<string, string> variables;
204 variables["[filename]"] = fnameRoot;
205 string pcoaDataFile = getOutputFileName("pcoa",variables);
206 util.openOutputFile(pcoaDataFile, pcaData);
207 pcaData.setf(ios::fixed, ios::floatfield);
208 pcaData.setf(ios::showpoint);
209 outputNames.push_back(pcoaDataFile);
210 outputTypes["pcoa"].push_back(pcoaDataFile);
211
212 ofstream pcaLoadings;
213 string loadingsFile = getOutputFileName("loadings",variables);
214 util.openOutputFile(loadingsFile, pcaLoadings);
215 pcaLoadings.setf(ios::fixed, ios::floatfield);
216 pcaLoadings.setf(ios::showpoint);
217 outputNames.push_back(loadingsFile);
218 outputTypes["loadings"].push_back(loadingsFile);
219
220 pcaLoadings << "axis\tloading\n";
221 for(int i=0;i<rank;i++){
222 pcaLoadings << i+1 << '\t' << d[i] * 100.0 / dsum << endl;
223 }
224
225 pcaData << "group";
226 for(int i=0;i<rank;i++){
227 pcaData << '\t' << "axis" << i+1;
228 }
229 pcaData << endl;
230
231 for(int i=0;i<rank;i++){
232 pcaData << name_list[i];
233 for(int j=0;j<rank;j++){
234 pcaData << '\t' << G[i][j];
235 }
236 pcaData << endl;
237 }
238 }
239 catch(exception& e) {
240 m->errorOut(e, "PCOACommand", "output");
241 exit(1);
242 }
243 }
244
245 /*********************************************************************************************************************************/
246
247