]> git.donarmstrong.com Git - mothur.git/blob - getgroupscommand.cpp
84b168a4c75164b339b7ed2b6f350d13021da5e4
[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->mothurOut(" contains " + toString(groupMap->getNumSeqs(Groups[i])) + " sequences."); m->mothurOutEndLine(); }
254                 m->mothurOutEndLine();
255                 
256                 
257                 
258                 if (outputNames.size() != 0) {
259                         m->mothurOutEndLine();
260                         m->mothurOut("Output File names: "); m->mothurOutEndLine();
261                         for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
262                         m->mothurOutEndLine();
263                         
264                         //set fasta file as new current fastafile
265                         string current = "";
266                         itTypes = outputTypes.find("fasta");
267                         if (itTypes != outputTypes.end()) {
268                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
269                         }
270                         
271                         itTypes = outputTypes.find("name");
272                         if (itTypes != outputTypes.end()) {
273                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setNameFile(current); }
274                         }
275                         
276                         itTypes = outputTypes.find("group");
277                         if (itTypes != outputTypes.end()) {
278                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setGroupFile(current); }
279                         }
280                         
281                         itTypes = outputTypes.find("list");
282                         if (itTypes != outputTypes.end()) {
283                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setListFile(current); }
284                         }
285                         
286                         itTypes = outputTypes.find("taxonomy");
287                         if (itTypes != outputTypes.end()) {
288                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTaxonomyFile(current); }
289                         }
290                 }
291                 
292                 return 0;               
293         }
294         
295         catch(exception& e) {
296                 m->errorOut(e, "GetGroupsCommand", "execute");
297                 exit(1);
298         }
299 }
300
301 //**********************************************************************************************************************
302 int GetGroupsCommand::readFasta(){
303         try {
304                 string thisOutputDir = outputDir;
305                 if (outputDir == "") {  thisOutputDir += m->hasPath(fastafile);  }
306                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(fastafile)) + "pick" + m->getExtension(fastafile);
307                 
308                 ofstream out;
309                 m->openOutputFile(outputFileName, out);
310                 
311                 ifstream in;
312                 m->openInputFile(fastafile, in);
313                 string name;
314                 
315                 bool wroteSomething = false;
316                 
317                 while(!in.eof()){
318                         if (m->control_pressed) { in.close();  out.close();  remove(outputFileName.c_str());  return 0; }
319                         
320                         Sequence currSeq(in);
321                         name = currSeq.getName();
322                         
323                         if (name != "") {
324                                 //if this name is in the accnos file
325                                 if (names.count(name) != 0) {
326                                         wroteSomething = true;
327                                         
328                                         currSeq.printSequence(out);
329                                 }
330                         }
331                         m->gobble(in);
332                 }
333                 in.close();     
334                 out.close();
335                 
336                 if (wroteSomething == false) {  m->mothurOut("Your file does NOT contain sequences from the groups you wish to get."); m->mothurOutEndLine();  }
337                 outputTypes["fasta"].push_back(outputFileName);  outputNames.push_back(outputFileName);
338                 
339                 return 0;
340                 
341         }
342         catch(exception& e) {
343                 m->errorOut(e, "GetGroupsCommand", "readFasta");
344                 exit(1);
345         }
346 }
347
348 //**********************************************************************************************************************
349 int GetGroupsCommand::readList(){
350         try {
351                 string thisOutputDir = outputDir;
352                 if (outputDir == "") {  thisOutputDir += m->hasPath(listfile);  }
353                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(listfile)) + "pick" +  m->getExtension(listfile);
354                 
355                 ofstream out;
356                 m->openOutputFile(outputFileName, out);
357                 
358                 ifstream in;
359                 m->openInputFile(listfile, in);
360                 
361                 bool wroteSomething = false;
362                 
363                 while(!in.eof()){
364                         //read in list vector
365                         ListVector list(in);
366                         
367                         //make a new list vector
368                         ListVector newList;
369                         newList.setLabel(list.getLabel());
370                         
371                         //for each bin
372                         for (int i = 0; i < list.getNumBins(); i++) {
373                                 if (m->control_pressed) { in.close();  out.close();  remove(outputFileName.c_str());  return 0; }
374                                 
375                                 //parse out names that are in accnos file
376                                 string binnames = list.get(i);
377                                 
378                                 string newNames = "";
379                                 while (binnames.find_first_of(',') != -1) { 
380                                         string name = binnames.substr(0,binnames.find_first_of(','));
381                                         binnames = binnames.substr(binnames.find_first_of(',')+1, binnames.length());
382                                         
383                                         //if that name is in the .accnos file, add it
384                                         if (names.count(name) != 0) {  newNames += name + ",";  }
385                                 }
386                                 
387                                 //get last name
388                                 if (names.count(binnames) != 0) {  newNames += binnames + ",";  }
389                                 
390                                 //if there are names in this bin add to new list
391                                 if (newNames != "") {  
392                                         newNames = newNames.substr(0, newNames.length()-1); //rip off extra comma
393                                         newList.push_back(newNames);    
394                                 }
395                         }
396                         
397                         //print new listvector
398                         if (newList.getNumBins() != 0) {
399                                 wroteSomething = true;
400                                 newList.print(out);
401                         }
402                         
403                         m->gobble(in);
404                 }
405                 in.close();     
406                 out.close();
407                 
408                 if (wroteSomething == false) {  m->mothurOut("Your file does NOT contain sequences from the groups you wish to get."); m->mothurOutEndLine();  }
409                 outputTypes["list"].push_back(outputFileName); outputNames.push_back(outputFileName);
410                 
411                 return 0;
412                 
413         }
414         catch(exception& e) {
415                 m->errorOut(e, "GetGroupsCommand", "readList");
416                 exit(1);
417         }
418 }
419 //**********************************************************************************************************************
420 int GetGroupsCommand::readName(){
421         try {
422                 string thisOutputDir = outputDir;
423                 if (outputDir == "") {  thisOutputDir += m->hasPath(namefile);  }
424                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(namefile)) + "pick" + m->getExtension(namefile);
425                 
426                 ofstream out;
427                 m->openOutputFile(outputFileName, out);
428                 
429                 ifstream in;
430                 m->openInputFile(namefile, in);
431                 string name, firstCol, secondCol;
432                 
433                 bool wroteSomething = false;
434                 
435                 while(!in.eof()){
436                         if (m->control_pressed) { in.close();  out.close();  remove(outputFileName.c_str());  return 0; }
437                         
438                         in >> firstCol;         m->gobble(in);          
439                         in >> secondCol;                        
440                         
441                         vector<string> parsedNames;
442                         m->splitAtComma(secondCol, parsedNames);
443                         
444                         vector<string> validSecond;  validSecond.clear();
445                         for (int i = 0; i < parsedNames.size(); i++) {
446                                 if (names.count(parsedNames[i]) != 0) {
447                                         validSecond.push_back(parsedNames[i]);
448                                 }
449                         }
450                         
451                         //if the name in the first column is in the set then print it and any other names in second column also in set
452                         if (names.count(firstCol) != 0) {
453                                 
454                                 wroteSomething = true;
455                                 
456                                 out << firstCol << '\t';
457                                 
458                                 //you know you have at least one valid second since first column is valid
459                                 for (int i = 0; i < validSecond.size()-1; i++) {  out << validSecond[i] << ',';  }
460                                 out << validSecond[validSecond.size()-1] << endl;
461                                 
462                                 //make first name in set you come to first column and then add the remaining names to second column
463                         }else {
464                                 
465                                 //you want part of this row
466                                 if (validSecond.size() != 0) {
467                                         
468                                         wroteSomething = true;
469                                         
470                                         out << validSecond[0] << '\t';
471                                         
472                                         //you know you have at least one valid second since first column is valid
473                                         for (int i = 0; i < validSecond.size()-1; i++) {  out << validSecond[i] << ',';  }
474                                         out << validSecond[validSecond.size()-1] << endl;
475                                 }
476                         }
477                         
478                         m->gobble(in);
479                 }
480                 in.close();
481                 out.close();
482                 
483                 if (wroteSomething == false) {  m->mothurOut("Your file does NOT contain sequences from the groups you wish to get."); m->mothurOutEndLine();  }
484                 outputTypes["name"].push_back(outputFileName); outputNames.push_back(outputFileName);
485                 
486                 return 0;
487         }
488         catch(exception& e) {
489                 m->errorOut(e, "GetGroupsCommand", "readName");
490                 exit(1);
491         }
492 }
493
494 //**********************************************************************************************************************
495 int GetGroupsCommand::readGroup(){
496         try {
497                 string thisOutputDir = outputDir;
498                 if (outputDir == "") {  thisOutputDir += m->hasPath(groupfile);  }
499                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(groupfile)) + "pick" + m->getExtension(groupfile);
500                 
501                 ofstream out;
502                 m->openOutputFile(outputFileName, out);
503                 
504                 ifstream in;
505                 m->openInputFile(groupfile, in);
506                 string name, group;
507                 
508                 bool wroteSomething = false;
509                 
510                 while(!in.eof()){
511                         if (m->control_pressed) { in.close();  out.close();  remove(outputFileName.c_str());  return 0; }
512                         
513                         in >> name;                             //read from first column
514                         in >> group;                    //read from second column
515                         
516                         //if this name is in the accnos file
517                         if (names.count(name) != 0) {
518                                 wroteSomething = true;
519                                 out << name << '\t' << group << endl;
520                         }
521                         
522                         m->gobble(in);
523                 }
524                 in.close();
525                 out.close();
526                 
527                 if (wroteSomething == false) {  m->mothurOut("Your file does NOT contain sequences from the groups you wish to get."); m->mothurOutEndLine();  }
528                 outputTypes["group"].push_back(outputFileName); outputNames.push_back(outputFileName);
529                 
530                 return 0;
531         }
532         catch(exception& e) {
533                 m->errorOut(e, "GetGroupsCommand", "readGroup");
534                 exit(1);
535         }
536 }
537 //**********************************************************************************************************************
538 int GetGroupsCommand::readTax(){
539         try {
540                 string thisOutputDir = outputDir;
541                 if (outputDir == "") {  thisOutputDir += m->hasPath(taxfile);  }
542                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(taxfile)) + "pick" + m->getExtension(taxfile);
543                 ofstream out;
544                 m->openOutputFile(outputFileName, out);
545                 
546                 ifstream in;
547                 m->openInputFile(taxfile, in);
548                 string name, tax;
549                 
550                 bool wroteSomething = false;
551                 
552                 while(!in.eof()){
553                         if (m->control_pressed) { in.close();  out.close();  remove(outputFileName.c_str());  return 0; }
554                         
555                         in >> name;                             //read from first column
556                         in >> tax;                      //read from second column
557                         
558                         //if this name is in the accnos file
559                         if (names.count(name) != 0) {
560                                 wroteSomething = true;
561                                 out << name << '\t' << tax << endl;
562                         }
563                         
564                         m->gobble(in);
565                 }
566                 in.close();
567                 out.close();
568                 
569                 if (wroteSomething == false) {  m->mothurOut("Your file does NOT contain sequences from the groups you wish to get."); m->mothurOutEndLine();  }
570                 outputTypes["taxonomy"].push_back(outputFileName); outputNames.push_back(outputFileName);
571                 
572                 return 0;
573         }
574         catch(exception& e) {
575                 m->errorOut(e, "GetGroupsCommand", "readTax");
576                 exit(1);
577         }
578 }
579 //**********************************************************************************************************************
580 void GetGroupsCommand::readAccnos(){
581         try {
582                 Groups.clear();
583                 
584                 ifstream in;
585                 m->openInputFile(accnosfile, in);
586                 string name;
587                 
588                 while(!in.eof()){
589                         in >> name;
590                         
591                         Groups.push_back(name);
592                         
593                         m->gobble(in);
594                 }
595                 in.close();             
596                 
597         }
598         catch(exception& e) {
599                 m->errorOut(e, "GetGroupsCommand", "readAccnos");
600                 exit(1);
601         }
602 }
603 //**********************************************************************************************************************
604 int GetGroupsCommand::fillNames(){
605         try {
606                 vector<string> seqs = groupMap->getNamesSeqs();
607                 
608                 for (int i = 0; i < seqs.size(); i++) {
609                         
610                         if (m->control_pressed) { return 0; }
611                         
612                         string group = groupMap->getGroup(seqs[i]);
613                         
614                         if (m->inUsersGroups(group, Groups)) {
615                                 names.insert(seqs[i]);
616                         }
617                 }
618                 
619                 return 0;
620         }
621         catch(exception& e) {
622                 m->errorOut(e, "GetGroupsCommand", "fillNames");
623                 exit(1);
624         }
625 }
626
627 //**********************************************************************************************************************
628
629