]> git.donarmstrong.com Git - mothur.git/blob - phylotypecommand.cpp
added citation function to commands
[mothur.git] / phylotypecommand.cpp
1 /*
2  *  phylotypecommand.cpp
3  *  Mothur
4  *
5  *  Created by westcott on 11/20/09.
6  *  Copyright 2009 Schloss Lab. All rights reserved.
7  *
8  */
9
10 #include "phylotypecommand.h"
11 #include "phylotree.h"
12 #include "listvector.hpp"
13 #include "rabundvector.hpp"
14 #include "sabundvector.hpp"
15
16 //**********************************************************************************************************************
17 vector<string> PhylotypeCommand::setParameters(){       
18         try {
19                 CommandParameter ptaxonomy("taxonomy", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(ptaxonomy);
20                 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
21                 CommandParameter pcutoff("cutoff", "Number", "", "-1", "", "", "",false,false); parameters.push_back(pcutoff);
22                 CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
23                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
24                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
25                 
26                 vector<string> myArray;
27                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
28                 return myArray;
29         }
30         catch(exception& e) {
31                 m->errorOut(e, "PhylotypeCommand", "setParameters");
32                 exit(1);
33         }
34 }
35 //**********************************************************************************************************************
36 string PhylotypeCommand::getHelpString(){       
37         try {
38                 string helpString = "";
39                 helpString += "The phylotype command reads a taxonomy file and outputs a .list, .rabund and .sabund file. \n";
40                 helpString += "The phylotype command parameter options are taxonomy, cutoff and label. The taxonomy parameter is required.\n";
41                 helpString += "The cutoff parameter allows you to specify the level you want to stop at.  The default is the highest level in your taxonomy file. \n";
42                 helpString += "For example: taxonomy = Bacteria;Bacteroidetes-Chlorobi;Bacteroidetes; - cutoff=2, would truncate the taxonomy to Bacteria;Bacteroidetes-Chlorobi; \n";
43                 helpString += "For the cutoff parameter levels count up from the root of the phylotree. This enables you to look at the grouping down to a specific resolution, say the genus level.\n";
44                 helpString += "The label parameter allows you to specify which level you would like, and are separated by dashes.  The default all levels in your taxonomy file. \n";
45                 helpString += "For the label parameter, levels count down from the root to keep the output similiar to mothur's other commands which report information from finer resolution to coarser resolutions.\n";
46                 helpString += "The phylotype command should be in the following format: \n";
47                 helpString += "phylotype(taxonomy=yourTaxonomyFile, cutoff=yourCutoff, label=yourLabels) \n";
48                 helpString += "Eaxample: phylotype(taxonomy=amazon.taxonomy, cutoff=5, label=1-3-5).\n";
49                 return helpString;
50         }
51         catch(exception& e) {
52                 m->errorOut(e, "PhylotypeCommand", "getHelpString");
53                 exit(1);
54         }
55 }
56
57 //**********************************************************************************************************************
58 PhylotypeCommand::PhylotypeCommand(){   
59         try {
60                 abort = true; calledHelp = true; 
61                 setParameters();
62                 vector<string> tempOutNames;
63                 outputTypes["list"] = tempOutNames;
64                 outputTypes["sabund"] = tempOutNames;
65                 outputTypes["rabund"] = tempOutNames;
66         }
67         catch(exception& e) {
68                 m->errorOut(e, "PhylotypeCommand", "PhylotypeCommand");
69                 exit(1);
70         }
71 }
72 /**********************************************************************************************************************/
73 PhylotypeCommand::PhylotypeCommand(string option)  {
74         try {
75                 abort = false; calledHelp = false;   
76                 
77                 //allow user to run help
78                 if(option == "help") { help(); abort = true; calledHelp = true; }
79                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
80                 
81                 else {
82                         vector<string> myArray = setParameters();
83                         
84                         OptionParser parser(option);
85                         map<string, string> parameters = parser.getParameters(); 
86                         
87                         ValidParameters validParameter;
88                         map<string, string>::iterator it;
89                         
90                         //check to make sure all parameters are valid for command
91                         for (it = parameters.begin(); it != parameters.end(); it++) { 
92                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
93                         }
94                         
95                         //initialize outputTypes
96                         vector<string> tempOutNames;
97                         outputTypes["list"] = tempOutNames;
98                         outputTypes["sabund"] = tempOutNames;
99                         outputTypes["rabund"] = tempOutNames;
100                         
101                         //if the user changes the input directory command factory will send this info to us in the output parameter 
102                         string inputDir = validParameter.validFile(parameters, "inputdir", false);              
103                         if (inputDir == "not found"){   inputDir = "";          }
104                         else {
105                                 string path;
106                                 it = parameters.find("taxonomy");
107                                 //user has given a template file
108                                 if(it != parameters.end()){ 
109                                         path = m->hasPath(it->second);
110                                         //if the user has not given a path then, add inputdir. else leave path alone.
111                                         if (path == "") {       parameters["taxonomy"] = inputDir + it->second;         }
112                                 }
113                                 
114                                 it = parameters.find("name");
115                                 //user has given a template file
116                                 if(it != parameters.end()){ 
117                                         path = m->hasPath(it->second);
118                                         //if the user has not given a path then, add inputdir. else leave path alone.
119                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
120                                 }
121                         }
122
123                         taxonomyFileName = validParameter.validFile(parameters, "taxonomy", true);
124                         if (taxonomyFileName == "not found") { 
125                                 taxonomyFileName = m->getTaxonomyFile(); 
126                                 if (taxonomyFileName != "") {  m->mothurOut("Using " + taxonomyFileName + " as input file for the taxonomy parameter."); m->mothurOutEndLine(); }
127                                 else { 
128                                         m->mothurOut("No valid current files. taxonomy is a required parameter."); m->mothurOutEndLine(); 
129                                         abort = true; 
130                                 }
131                         }else if (taxonomyFileName == "not open") { abort = true; }     
132                         
133                         namefile = validParameter.validFile(parameters, "name", true);
134                         if (namefile == "not open") { abort = true; }
135                         else if (namefile == "not found") { namefile = ""; }
136                         else { readNamesFile(); }       
137                         
138                         //if the user changes the output directory command factory will send this info to us in the output parameter 
139                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  
140                                 outputDir = ""; 
141                                 outputDir += m->hasPath(taxonomyFileName); //if user entered a file with a path then preserve it        
142                         }
143                         
144                         string temp = validParameter.validFile(parameters, "cutoff", false);
145                         if (temp == "not found") { temp = "-1"; }
146                         convert(temp, cutoff); 
147                         
148                         label = validParameter.validFile(parameters, "label", false);                   
149                         if (label == "not found") { label = ""; allLines = 1; }
150                         else { 
151                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
152                                 else { allLines = 1;  }
153                         }
154                         
155                 }
156         }
157         catch(exception& e) {
158                 m->errorOut(e, "PhylotypeCommand", "PhylotypeCommand");
159                 exit(1);
160         }
161 }
162 /**********************************************************************************************************************/
163
164 int PhylotypeCommand::execute(){
165         try {
166         
167                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
168                 
169                 //reads in taxonomy file and makes all the taxonomies the same length 
170                 //by appending the last taxon to a given taxonomy as many times as needed to 
171                 //make it as long as the longest taxonomy in the file 
172                 TaxEqualizer* taxEqual = new TaxEqualizer(taxonomyFileName, cutoff, outputDir);
173                 
174                 if (m->control_pressed) { delete taxEqual; return 0; }
175                 
176                 string equalizedTaxFile = taxEqual->getEqualizedTaxFile();
177                 
178                 delete taxEqual;
179                 
180                 //build taxonomy tree from equalized file
181                 PhyloTree* tree = new PhyloTree(equalizedTaxFile);
182                 vector<int> leaves = tree->getGenusNodes();
183                 
184                 //store leaf nodes in current map
185                 for (int i = 0; i < leaves.size(); i++)         {       currentNodes[leaves[i]] = leaves[i];    }
186                 
187                 bool done = false;
188                 if (tree->get(leaves[0]).parent == -1) {  m->mothurOut("Empty Tree"); m->mothurOutEndLine();    done = true;    }
189                 
190                 if (m->control_pressed) { delete tree; return 0; }
191                 
192                 string fileroot = outputDir + m->getRootName(m->getSimpleName(taxonomyFileName));
193                 
194                 ofstream outList;
195                 string outputListFile = fileroot + "tx.list";
196                 m->openOutputFile(outputListFile, outList);
197                 ofstream outSabund;
198                 string outputSabundFile = fileroot + "tx.sabund";
199                 m->openOutputFile(outputSabundFile, outSabund);
200                 ofstream outRabund;
201                 string outputRabundFile = fileroot + "tx.rabund";
202                 m->openOutputFile(outputRabundFile, outRabund);
203                 
204                 outputNames.push_back(outputListFile); outputTypes["list"].push_back(outputListFile);
205                 outputNames.push_back(outputSabundFile); outputTypes["sabund"].push_back(outputSabundFile);
206                 outputNames.push_back(outputRabundFile); outputTypes["rabund"].push_back(outputRabundFile);
207                 
208                 int count = 1;          
209                 //start at leaves of tree and work towards root, processing the labels the user wants
210                 while((!done) && ((allLines == 1) || (labels.size() != 0))) {
211                 
212                         string level = toString(count); 
213                         count++;
214                         
215                         if (m->control_pressed) { 
216                                 outRabund.close(); outSabund.close(); outList.close();
217                                 for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  }
218                                 delete tree; return 0; 
219                         }
220                         
221                         //is this a level the user want output for
222                         if(allLines == 1 || labels.count(level) == 1){  
223                                 
224                                 //output level
225                                 m->mothurOut(level); m->mothurOutEndLine();
226                                 
227                                 ListVector list;
228                                 list.setLabel(level);
229                                 //go through nodes and build listvector 
230                                 for (itCurrent = currentNodes.begin(); itCurrent != currentNodes.end(); itCurrent++) {
231                         
232                                         //get parents
233                                         TaxNode node = tree->get(itCurrent->first);
234                                         parentNodes[node.parent] = node.parent;
235                                         
236                                         vector<string> names = node.accessions;
237                                         
238                                         //make the names compatable with listvector
239                                         string name = "";
240                                         for (int i = 0; i < names.size(); i++) {  
241                                                 if (namefile != "") {   
242                                                         map<string, string>::iterator itNames = namemap.find(names[i]);  //make sure this name is in namefile
243                 
244                                                         if (itNames != namemap.end()) {  name += namemap[names[i]] + ",";   } //you found it in namefile
245                                                         else { m->mothurOut(names[i] + " is not in your namefile, please correct."); m->mothurOutEndLine(); exit(1);  }
246                                                         
247                                                 }else{   name += names[i] + ",";        }
248                                         }
249                                         name = name.substr(0, name.length()-1);  //rip off extra ','
250                                         
251                                         //add bin to list vector
252                                         list.push_back(name);
253                                 }       
254                                 
255                                 //print listvector
256                                 list.print(outList);
257                                 //print rabund
258                                 list.getRAbundVector().print(outRabund);
259                                 //print sabund
260                                 list.getSAbundVector().print(outSabund);
261                         
262                                 labels.erase(level);
263                                 
264                         }else {
265                                 
266                                 //just get parents
267                                 for (itCurrent = currentNodes.begin(); itCurrent != currentNodes.end(); itCurrent++) {
268                                         int parent = tree->get(itCurrent->first).parent;
269                                         parentNodes[parent] = parent;
270                                 }
271                         }
272                         
273                         //move up a level
274                         currentNodes = parentNodes;
275                         parentNodes.clear();
276                         
277                         //have we reached the rootnode
278                         if (tree->get(currentNodes.begin()->first).parent == -1)  {  done = true;  }
279                 }
280                         
281                 outList.close();
282                 outSabund.close();
283                 outRabund.close();      
284                 
285                 delete tree;
286                 
287                 if (m->control_pressed) { 
288                         for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  }
289                         return 0; 
290                 }
291                 
292                 //set list file as new current listfile
293                 string current = "";
294                 itTypes = outputTypes.find("list");
295                 if (itTypes != outputTypes.end()) {
296                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setListFile(current); }
297                 }
298                 
299                 //set rabund file as new current rabundfile
300                 itTypes = outputTypes.find("rabund");
301                 if (itTypes != outputTypes.end()) {
302                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setRabundFile(current); }
303                 }
304                 
305                 //set sabund file as new current sabundfile
306                 itTypes = outputTypes.find("sabund");
307                 if (itTypes != outputTypes.end()) {
308                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setSabundFile(current); }
309                 }               
310                 
311                 m->mothurOutEndLine();
312                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
313                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
314                 m->mothurOutEndLine();
315                 
316                 return 0;               
317         }
318
319         catch(exception& e) {
320                 m->errorOut(e, "PhylotypeCommand", "execute");
321                 exit(1);
322         }
323 }
324 /*****************************************************************/
325 int PhylotypeCommand::readNamesFile() {
326         try {
327                                 
328                 ifstream in;
329                 m->openInputFile(namefile, in);
330                 
331                 string first, second;
332                 map<string, string>::iterator itNames;
333                 
334                 while(!in.eof()) {
335                         in >> first >> second; m->gobble(in);
336                         
337                         itNames = namemap.find(first);
338                         if (itNames == namemap.end()) {  
339                                 namemap[first] = second; 
340                         }else {  m->mothurOut(first + " has already been seen in namefile, disregarding names file."); m->mothurOutEndLine(); in.close(); namemap.clear(); namefile = ""; return 1; }                   
341                 }
342                 in.close();
343                 
344                 return 0;
345         }
346         catch(exception& e) {
347                 m->errorOut(e, "PhylotypeCommand", "readNamesFile");
348                 exit(1);
349         }
350 }
351
352 /**********************************************************************************************************************/