]> git.donarmstrong.com Git - mothur.git/blob - getgroupscommand.cpp
added set.current and get.current commands and modified existing commands to update...
[mothur.git] / getgroupscommand.cpp
1 /*
2  *  getgroupscommand.cpp
3  *  Mothur
4  *
5  *  Created by westcott on 11/10/10.
6  *  Copyright 2010 Schloss Lab. All rights reserved.
7  *
8  */
9
10 #include "getgroupscommand.h"
11 #include "sequence.hpp"
12 #include "listvector.hpp"
13 #include "sharedutilities.h"
14
15 //**********************************************************************************************************************
16 vector<string> GetGroupsCommand::getValidParameters(){  
17         try {
18                 string Array[] =  {"fasta","name", "group", "accnos", "groups","list","taxonomy","outputdir","inputdir" };
19                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
20                 return myArray;
21         }
22         catch(exception& e) {
23                 m->errorOut(e, "GetGroupsCommand", "getValidParameters");
24                 exit(1);
25         }
26 }
27 //**********************************************************************************************************************
28 GetGroupsCommand::GetGroupsCommand(){   
29         try {
30                 abort = true; calledHelp = true; 
31                 vector<string> tempOutNames;
32                 outputTypes["fasta"] = tempOutNames;
33                 outputTypes["taxonomy"] = tempOutNames;
34                 outputTypes["name"] = tempOutNames;
35                 outputTypes["group"] = tempOutNames;
36                 outputTypes["list"] = tempOutNames;
37         }
38         catch(exception& e) {
39                 m->errorOut(e, "GetGroupsCommand", "GetGroupsCommand");
40                 exit(1);
41         }
42 }
43 //**********************************************************************************************************************
44 vector<string> GetGroupsCommand::getRequiredParameters(){       
45         try {
46                 string Array[] =  {"group"};
47                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
48                 return myArray;
49         }
50         catch(exception& e) {
51                 m->errorOut(e, "GetGroupsCommand", "getRequiredParameters");
52                 exit(1);
53         }
54 }
55 //**********************************************************************************************************************
56 vector<string> GetGroupsCommand::getRequiredFiles(){    
57         try {
58                 vector<string> myArray;
59                 return myArray;
60         }
61         catch(exception& e) {
62                 m->errorOut(e, "GetGroupsCommand", "getRequiredFiles");
63                 exit(1);
64         }
65 }
66 //**********************************************************************************************************************
67 GetGroupsCommand::GetGroupsCommand(string option)  {
68         try {
69                 abort = false; calledHelp = false;   
70                 
71                 //allow user to run help
72                 if(option == "help") { help(); abort = true; calledHelp = true; }
73                 
74                 else {
75                         //valid paramters for this command
76                         string Array[] =  {"fasta","name", "group", "accnos", "groups", "list","taxonomy","outputdir","inputdir" };
77                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
78                         
79                         OptionParser parser(option);
80                         map<string,string> parameters = parser.getParameters();
81                         
82                         ValidParameters validParameter;
83                         map<string,string>::iterator it;
84                         
85                         //check to make sure all parameters are valid for command
86                         for (it = parameters.begin(); it != parameters.end(); it++) { 
87                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
88                         }
89                         
90                         //initialize outputTypes
91                         vector<string> tempOutNames;
92                         outputTypes["fasta"] = tempOutNames;
93                         outputTypes["taxonomy"] = tempOutNames;
94                         outputTypes["name"] = tempOutNames;
95                         outputTypes["group"] = tempOutNames;
96                         outputTypes["list"] = tempOutNames;
97                         
98                         
99                         //if the user changes the output directory command factory will send this info to us in the output parameter 
100                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = "";         }
101                         
102                         //if the user changes the input directory command factory will send this info to us in the output parameter 
103                         string inputDir = validParameter.validFile(parameters, "inputdir", false);              
104                         if (inputDir == "not found"){   inputDir = "";          }
105                         else {
106                                 string path;
107                                 it = parameters.find("fasta");
108                                 //user has given a template file
109                                 if(it != parameters.end()){ 
110                                         path = m->hasPath(it->second);
111                                         //if the user has not given a path then, add inputdir. else leave path alone.
112                                         if (path == "") {       parameters["fasta"] = inputDir + it->second;            }
113                                 }
114                                 
115                                 it = parameters.find("accnos");
116                                 //user has given a template file
117                                 if(it != parameters.end()){ 
118                                         path = m->hasPath(it->second);
119                                         //if the user has not given a path then, add inputdir. else leave path alone.
120                                         if (path == "") {       parameters["accnos"] = inputDir + it->second;           }
121                                 }
122                                 
123                                 it = parameters.find("list");
124                                 //user has given a template file
125                                 if(it != parameters.end()){ 
126                                         path = m->hasPath(it->second);
127                                         //if the user has not given a path then, add inputdir. else leave path alone.
128                                         if (path == "") {       parameters["list"] = inputDir + it->second;             }
129                                 }
130                                 
131                                 it = parameters.find("name");
132                                 //user has given a template file
133                                 if(it != parameters.end()){ 
134                                         path = m->hasPath(it->second);
135                                         //if the user has not given a path then, add inputdir. else leave path alone.
136                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
137                                 }
138                                 
139                                 it = parameters.find("group");
140                                 //user has given a template file
141                                 if(it != parameters.end()){ 
142                                         path = m->hasPath(it->second);
143                                         //if the user has not given a path then, add inputdir. else leave path alone.
144                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
145                                 }
146                                 
147                                 it = parameters.find("taxonomy");
148                                 //user has given a template file
149                                 if(it != parameters.end()){ 
150                                         path = m->hasPath(it->second);
151                                         //if the user has not given a path then, add inputdir. else leave path alone.
152                                         if (path == "") {       parameters["taxonomy"] = inputDir + it->second;         }
153                                 }
154                         }
155                         
156                         
157                         //check for required parameters
158                         accnosfile = validParameter.validFile(parameters, "accnos", true);
159                         if (accnosfile == "not open") { abort = true; }
160                         else if (accnosfile == "not found") {  accnosfile = ""; }       
161                         
162                         fastafile = validParameter.validFile(parameters, "fasta", true);
163                         if (fastafile == "not open") { abort = true; }
164                         else if (fastafile == "not found") {  fastafile = "";  }        
165                         
166                         namefile = validParameter.validFile(parameters, "name", true);
167                         if (namefile == "not open") { abort = true; }
168                         else if (namefile == "not found") {  namefile = "";  }  
169                         
170                         groupfile = validParameter.validFile(parameters, "group", true);
171                         if (groupfile == "not open") { abort = true; }
172                         else if (groupfile == "not found") {  groupfile = "";  m->mothurOut("You must provide a group file."); m->mothurOutEndLine(); abort = true; }   
173                         
174                         listfile = validParameter.validFile(parameters, "list", true);
175                         if (listfile == "not open") { abort = true; }
176                         else if (listfile == "not found") {  listfile = "";  }
177                         
178                         taxfile = validParameter.validFile(parameters, "taxonomy", true);
179                         if (taxfile == "not open") { abort = true; }
180                         else if (taxfile == "not found") {  taxfile = "";  }
181                         
182                         groups = validParameter.validFile(parameters, "groups", false);                 
183                         if (groups == "not found") { groups = ""; }
184                         else { 
185                                 m->splitAtDash(groups, Groups);
186                         }
187                         
188                         if ((accnosfile == "") && (Groups.size() == 0)) { m->mothurOut("You must provide an accnos file or specify groups using the groups parameter."); m->mothurOutEndLine(); abort = true; }
189                         
190                         if ((fastafile == "") && (namefile == "") && (groupfile == "")  && (listfile == "") && (taxfile == ""))  { m->mothurOut("You must provide at least one of the following: fasta, name, taxonomy, group or list."); m->mothurOutEndLine(); abort = true; }
191                 }
192                 
193         }
194         catch(exception& e) {
195                 m->errorOut(e, "GetGroupsCommand", "GetGroupsCommand");
196                 exit(1);
197         }
198 }
199 //**********************************************************************************************************************
200
201 void GetGroupsCommand::help(){
202         try {
203                 m->mothurOut("The get.groups command selects sequences from a specfic group or set of groups from the following file types: fasta, name, group, list, taxonomy.\n");
204                 m->mothurOut("It outputs a file containing the sequences in the those specified groups.\n");
205                 m->mothurOut("The get.groups command parameters are accnos, fasta, name, group, list, taxonomy and groups. The group parameter is required.\n");
206                 m->mothurOut("You must also provide an accnos containing the list of groups to get or set the groups parameter to the groups you wish to select.\n");
207                 m->mothurOut("The groups parameter allows you to specify which of the groups in your groupfile you would like.  You can separate group names with dashes.\n");
208                 m->mothurOut("The get.groups command should be in the following format: get.groups(accnos=yourAccnos, fasta=yourFasta, group=yourGroupFile).\n");
209                 m->mothurOut("Example get.groups(accnos=amazon.accnos, fasta=amazon.fasta, group=amazon.groups).\n");
210                 m->mothurOut("or get.groups(groups=pasture, fasta=amazon.fasta, group=amazon.groups).\n");
211                 m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n");
212         }
213         catch(exception& e) {
214                 m->errorOut(e, "GetGroupsCommand", "help");
215                 exit(1);
216         }
217 }
218
219 //**********************************************************************************************************************
220
221 int GetGroupsCommand::execute(){
222         try {
223                 
224                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
225                 
226                 groupMap = new GroupMap(groupfile);
227                 groupMap->readMap();
228                 
229                 //get groups you want to remove
230                 if (accnosfile != "") { readAccnos(); }
231                 
232                 //make sure groups are valid
233                 //takes care of user setting groupNames that are invalid or setting groups=all
234                 SharedUtil* util = new SharedUtil();
235                 util->setGroups(Groups, groupMap->namesOfGroups);
236                 delete util;
237                 
238                 //fill names with names of sequences that are from the groups we want to remove 
239                 fillNames();
240                 
241                 if (m->control_pressed) { delete groupMap; return 0; }
242                 
243                 //read through the correct file and output lines you want to keep
244                 if (namefile != "")                     {               readName();             }
245                 if (fastafile != "")            {               readFasta();    }
246                 if (groupfile != "")            {               readGroup();    }
247                 if (listfile != "")                     {               readList();             }
248                 if (taxfile != "")                      {               readTax();              }
249                 
250                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str()); } return 0; }
251                 
252                 m->mothurOut("Selected " + toString(names.size()) + " sequences. From the groups: "); m->mothurOutEndLine();
253                 for (int i = 0; i < Groups.size(); i++) {       m->mothurOut(Groups[i]); m->mothurOutEndLine(); }
254                 m->mothurOutEndLine();
255                 
256                 if (outputNames.size() != 0) {
257                         m->mothurOutEndLine();
258                         m->mothurOut("Output File names: "); m->mothurOutEndLine();
259                         for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
260                         m->mothurOutEndLine();
261                         
262                         //set fasta file as new current fastafile
263                         string current = "";
264                         itTypes = outputTypes.find("fasta");
265                         if (itTypes != outputTypes.end()) {
266                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
267                         }
268                         
269                         itTypes = outputTypes.find("name");
270                         if (itTypes != outputTypes.end()) {
271                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setNameFile(current); }
272                         }
273                         
274                         itTypes = outputTypes.find("group");
275                         if (itTypes != outputTypes.end()) {
276                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setGroupFile(current); }
277                         }
278                         
279                         itTypes = outputTypes.find("list");
280                         if (itTypes != outputTypes.end()) {
281                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setListFile(current); }
282                         }
283                         
284                         itTypes = outputTypes.find("taxonomy");
285                         if (itTypes != outputTypes.end()) {
286                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTaxonomyFile(current); }
287                         }
288                 }
289                 
290                 return 0;               
291         }
292         
293         catch(exception& e) {
294                 m->errorOut(e, "GetGroupsCommand", "execute");
295                 exit(1);
296         }
297 }
298
299 //**********************************************************************************************************************
300 int GetGroupsCommand::readFasta(){
301         try {
302                 string thisOutputDir = outputDir;
303                 if (outputDir == "") {  thisOutputDir += m->hasPath(fastafile);  }
304                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(fastafile)) + "pick" + m->getExtension(fastafile);
305                 
306                 ofstream out;
307                 m->openOutputFile(outputFileName, out);
308                 
309                 ifstream in;
310                 m->openInputFile(fastafile, in);
311                 string name;
312                 
313                 bool wroteSomething = false;
314                 
315                 while(!in.eof()){
316                         if (m->control_pressed) { in.close();  out.close();  remove(outputFileName.c_str());  return 0; }
317                         
318                         Sequence currSeq(in);
319                         name = currSeq.getName();
320                         
321                         if (name != "") {
322                                 //if this name is in the accnos file
323                                 if (names.count(name) != 0) {
324                                         wroteSomething = true;
325                                         
326                                         currSeq.printSequence(out);
327                                 }
328                         }
329                         m->gobble(in);
330                 }
331                 in.close();     
332                 out.close();
333                 
334                 if (wroteSomething == false) {  m->mothurOut("Your file does NOT contain sequences from the groups you wish to get."); m->mothurOutEndLine();  }
335                 outputTypes["fasta"].push_back(outputFileName);  outputNames.push_back(outputFileName);
336                 
337                 return 0;
338                 
339         }
340         catch(exception& e) {
341                 m->errorOut(e, "GetGroupsCommand", "readFasta");
342                 exit(1);
343         }
344 }
345
346 //**********************************************************************************************************************
347 int GetGroupsCommand::readList(){
348         try {
349                 string thisOutputDir = outputDir;
350                 if (outputDir == "") {  thisOutputDir += m->hasPath(listfile);  }
351                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(listfile)) + "pick" +  m->getExtension(listfile);
352                 
353                 ofstream out;
354                 m->openOutputFile(outputFileName, out);
355                 
356                 ifstream in;
357                 m->openInputFile(listfile, in);
358                 
359                 bool wroteSomething = false;
360                 
361                 while(!in.eof()){
362                         //read in list vector
363                         ListVector list(in);
364                         
365                         //make a new list vector
366                         ListVector newList;
367                         newList.setLabel(list.getLabel());
368                         
369                         //for each bin
370                         for (int i = 0; i < list.getNumBins(); i++) {
371                                 if (m->control_pressed) { in.close();  out.close();  remove(outputFileName.c_str());  return 0; }
372                                 
373                                 //parse out names that are in accnos file
374                                 string binnames = list.get(i);
375                                 
376                                 string newNames = "";
377                                 while (binnames.find_first_of(',') != -1) { 
378                                         string name = binnames.substr(0,binnames.find_first_of(','));
379                                         binnames = binnames.substr(binnames.find_first_of(',')+1, binnames.length());
380                                         
381                                         //if that name is in the .accnos file, add it
382                                         if (names.count(name) != 0) {  newNames += name + ",";  }
383                                 }
384                                 
385                                 //get last name
386                                 if (names.count(binnames) != 0) {  newNames += binnames + ",";  }
387                                 
388                                 //if there are names in this bin add to new list
389                                 if (newNames != "") {  
390                                         newNames = newNames.substr(0, newNames.length()-1); //rip off extra comma
391                                         newList.push_back(newNames);    
392                                 }
393                         }
394                         
395                         //print new listvector
396                         if (newList.getNumBins() != 0) {
397                                 wroteSomething = true;
398                                 newList.print(out);
399                         }
400                         
401                         m->gobble(in);
402                 }
403                 in.close();     
404                 out.close();
405                 
406                 if (wroteSomething == false) {  m->mothurOut("Your file does NOT contain sequences from the groups you wish to get."); m->mothurOutEndLine();  }
407                 outputTypes["list"].push_back(outputFileName); outputNames.push_back(outputFileName);
408                 
409                 return 0;
410                 
411         }
412         catch(exception& e) {
413                 m->errorOut(e, "GetGroupsCommand", "readList");
414                 exit(1);
415         }
416 }
417 //**********************************************************************************************************************
418 int GetGroupsCommand::readName(){
419         try {
420                 string thisOutputDir = outputDir;
421                 if (outputDir == "") {  thisOutputDir += m->hasPath(namefile);  }
422                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(namefile)) + "pick" + m->getExtension(namefile);
423                 
424                 ofstream out;
425                 m->openOutputFile(outputFileName, out);
426                 
427                 ifstream in;
428                 m->openInputFile(namefile, in);
429                 string name, firstCol, secondCol;
430                 
431                 bool wroteSomething = false;
432                 
433                 while(!in.eof()){
434                         if (m->control_pressed) { in.close();  out.close();  remove(outputFileName.c_str());  return 0; }
435                         
436                         in >> firstCol;         m->gobble(in);          
437                         in >> secondCol;                        
438                         
439                         vector<string> parsedNames;
440                         m->splitAtComma(secondCol, parsedNames);
441                         
442                         vector<string> validSecond;  validSecond.clear();
443                         for (int i = 0; i < parsedNames.size(); i++) {
444                                 if (names.count(parsedNames[i]) != 0) {
445                                         validSecond.push_back(parsedNames[i]);
446                                 }
447                         }
448                         
449                         //if the name in the first column is in the set then print it and any other names in second column also in set
450                         if (names.count(firstCol) != 0) {
451                                 
452                                 wroteSomething = true;
453                                 
454                                 out << firstCol << '\t';
455                                 
456                                 //you know you have at least one valid second since first column is valid
457                                 for (int i = 0; i < validSecond.size()-1; i++) {  out << validSecond[i] << ',';  }
458                                 out << validSecond[validSecond.size()-1] << endl;
459                                 
460                                 //make first name in set you come to first column and then add the remaining names to second column
461                         }else {
462                                 
463                                 //you want part of this row
464                                 if (validSecond.size() != 0) {
465                                         
466                                         wroteSomething = true;
467                                         
468                                         out << validSecond[0] << '\t';
469                                         
470                                         //you know you have at least one valid second since first column is valid
471                                         for (int i = 0; i < validSecond.size()-1; i++) {  out << validSecond[i] << ',';  }
472                                         out << validSecond[validSecond.size()-1] << endl;
473                                 }
474                         }
475                         
476                         m->gobble(in);
477                 }
478                 in.close();
479                 out.close();
480                 
481                 if (wroteSomething == false) {  m->mothurOut("Your file does NOT contain sequences from the groups you wish to get."); m->mothurOutEndLine();  }
482                 outputTypes["name"].push_back(outputFileName); outputNames.push_back(outputFileName);
483                 
484                 return 0;
485         }
486         catch(exception& e) {
487                 m->errorOut(e, "GetGroupsCommand", "readName");
488                 exit(1);
489         }
490 }
491
492 //**********************************************************************************************************************
493 int GetGroupsCommand::readGroup(){
494         try {
495                 string thisOutputDir = outputDir;
496                 if (outputDir == "") {  thisOutputDir += m->hasPath(groupfile);  }
497                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(groupfile)) + "pick" + m->getExtension(groupfile);
498                 
499                 ofstream out;
500                 m->openOutputFile(outputFileName, out);
501                 
502                 ifstream in;
503                 m->openInputFile(groupfile, in);
504                 string name, group;
505                 
506                 bool wroteSomething = false;
507                 
508                 while(!in.eof()){
509                         if (m->control_pressed) { in.close();  out.close();  remove(outputFileName.c_str());  return 0; }
510                         
511                         in >> name;                             //read from first column
512                         in >> group;                    //read from second column
513                         
514                         //if this name is in the accnos file
515                         if (names.count(name) != 0) {
516                                 wroteSomething = true;
517                                 out << name << '\t' << group << endl;
518                         }
519                         
520                         m->gobble(in);
521                 }
522                 in.close();
523                 out.close();
524                 
525                 if (wroteSomething == false) {  m->mothurOut("Your file does NOT contain sequences from the groups you wish to get."); m->mothurOutEndLine();  }
526                 outputTypes["group"].push_back(outputFileName); outputNames.push_back(outputFileName);
527                 
528                 return 0;
529         }
530         catch(exception& e) {
531                 m->errorOut(e, "GetGroupsCommand", "readGroup");
532                 exit(1);
533         }
534 }
535 //**********************************************************************************************************************
536 int GetGroupsCommand::readTax(){
537         try {
538                 string thisOutputDir = outputDir;
539                 if (outputDir == "") {  thisOutputDir += m->hasPath(taxfile);  }
540                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(taxfile)) + "pick" + m->getExtension(taxfile);
541                 ofstream out;
542                 m->openOutputFile(outputFileName, out);
543                 
544                 ifstream in;
545                 m->openInputFile(taxfile, in);
546                 string name, tax;
547                 
548                 bool wroteSomething = false;
549                 
550                 while(!in.eof()){
551                         if (m->control_pressed) { in.close();  out.close();  remove(outputFileName.c_str());  return 0; }
552                         
553                         in >> name;                             //read from first column
554                         in >> tax;                      //read from second column
555                         
556                         //if this name is in the accnos file
557                         if (names.count(name) != 0) {
558                                 wroteSomething = true;
559                                 out << name << '\t' << tax << endl;
560                         }
561                         
562                         m->gobble(in);
563                 }
564                 in.close();
565                 out.close();
566                 
567                 if (wroteSomething == false) {  m->mothurOut("Your file does NOT contain sequences from the groups you wish to get."); m->mothurOutEndLine();  }
568                 outputTypes["taxonomy"].push_back(outputFileName); outputNames.push_back(outputFileName);
569                 
570                 return 0;
571         }
572         catch(exception& e) {
573                 m->errorOut(e, "GetGroupsCommand", "readTax");
574                 exit(1);
575         }
576 }
577 //**********************************************************************************************************************
578 void GetGroupsCommand::readAccnos(){
579         try {
580                 Groups.clear();
581                 
582                 ifstream in;
583                 m->openInputFile(accnosfile, in);
584                 string name;
585                 
586                 while(!in.eof()){
587                         in >> name;
588                         
589                         Groups.push_back(name);
590                         
591                         m->gobble(in);
592                 }
593                 in.close();             
594                 
595         }
596         catch(exception& e) {
597                 m->errorOut(e, "GetGroupsCommand", "readAccnos");
598                 exit(1);
599         }
600 }
601 //**********************************************************************************************************************
602 int GetGroupsCommand::fillNames(){
603         try {
604                 vector<string> seqs = groupMap->getNamesSeqs();
605                 
606                 for (int i = 0; i < seqs.size(); i++) {
607                         
608                         if (m->control_pressed) { return 0; }
609                         
610                         string group = groupMap->getGroup(seqs[i]);
611                         
612                         if (m->inUsersGroups(group, Groups)) {
613                                 names.insert(seqs[i]);
614                         }
615                 }
616                 
617                 return 0;
618         }
619         catch(exception& e) {
620                 m->errorOut(e, "GetGroupsCommand", "fillNames");
621                 exit(1);
622         }
623 }
624
625 //**********************************************************************************************************************
626
627