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