]> git.donarmstrong.com Git - mothur.git/blob - listseqscommand.cpp
working on pam
[mothur.git] / listseqscommand.cpp
1 /*
2  *  listseqscommand.cpp
3  *  Mothur
4  *
5  *  Created by Sarah Westcott on 7/8/09.
6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
7  *
8  */
9
10 #include "listseqscommand.h"
11 #include "sequence.hpp"
12 #include "listvector.hpp"
13 #include "counttable.h"
14
15
16 //**********************************************************************************************************************
17 vector<string> ListSeqsCommand::setParameters(){        
18         try {
19         CommandParameter pfastq("fastq", "InputTypes", "", "", "FNGLT", "FNGLT", "none","accnos",false,false,true); parameters.push_back(pfastq);
20                 CommandParameter pfasta("fasta", "InputTypes", "", "", "FNGLT", "FNGLT", "none","accnos",false,false,true); parameters.push_back(pfasta);
21                 CommandParameter pname("name", "InputTypes", "", "", "FNGLT", "FNGLT", "none","accnos",false,false,true); parameters.push_back(pname);
22         CommandParameter pcount("count", "InputTypes", "", "", "FNGLT", "FNGLT", "none","accnos",false,false,true); parameters.push_back(pcount);
23                 CommandParameter pgroup("group", "InputTypes", "", "", "FNGLT", "FNGLT", "none","accnos",false,false,true); parameters.push_back(pgroup);
24                 CommandParameter plist("list", "InputTypes", "", "", "FNGLT", "FNGLT", "none","accnos",false,false,true); parameters.push_back(plist);
25                 CommandParameter ptaxonomy("taxonomy", "InputTypes", "", "", "FNGLT", "FNGLT", "none","accnos",false,false,true); parameters.push_back(ptaxonomy);
26                 CommandParameter palignreport("alignreport", "InputTypes", "", "", "FNGLT", "FNGLT", "none","accnos",false,false); parameters.push_back(palignreport);
27                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
28                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
29                 
30                 vector<string> myArray;
31                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
32                 return myArray;
33         }
34         catch(exception& e) {
35                 m->errorOut(e, "ListSeqsCommand", "setParameters");
36                 exit(1);
37         }
38 }
39 //**********************************************************************************************************************
40 string ListSeqsCommand::getHelpString(){        
41         try {
42                 string helpString = "";
43                 helpString += "The list.seqs command reads a fasta, name, group, count, list, taxonomy, fastq or alignreport file and outputs a .accnos file containing sequence names.\n";
44                 helpString += "The list.seqs command parameters are fasta, name, group, count, list, taxonomy, fastq and alignreport.  You must provide one of these parameters.\n";
45                 helpString += "The list.seqs command should be in the following format: list.seqs(fasta=yourFasta).\n";
46                 helpString += "Example list.seqs(fasta=amazon.fasta).\n";
47                 helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n";
48                 return helpString;
49         }
50         catch(exception& e) {
51                 m->errorOut(e, "ListSeqsCommand", "getHelpString");
52                 exit(1);
53         }
54 }
55 //**********************************************************************************************************************
56 string ListSeqsCommand::getOutputPattern(string type) {
57     try {
58         string pattern = "";
59         
60         if (type == "accnos") {  pattern = "[filename],accnos"; } 
61         else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
62         
63         return pattern;
64     }
65     catch(exception& e) {
66         m->errorOut(e, "ListSeqsCommand", "getOutputPattern");
67         exit(1);
68     }
69 }
70 //**********************************************************************************************************************
71 ListSeqsCommand::ListSeqsCommand(){     
72         try {
73                 abort = true; calledHelp = true; 
74                 setParameters();
75                 vector<string> tempOutNames;
76                 outputTypes["accnos"] = tempOutNames;
77         }
78         catch(exception& e) {
79                 m->errorOut(e, "ListSeqsCommand", "ListSeqsCommand");
80                 exit(1);
81         }
82 }
83 //**********************************************************************************************************************
84
85 ListSeqsCommand::ListSeqsCommand(string option)  {
86         try {
87                 abort = false; calledHelp = false;   
88                 
89                 //allow user to run help
90                 if(option == "help") { help(); abort = true; calledHelp = true; }
91                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
92                 else {
93                         vector<string> myArray = setParameters();
94                         
95                         OptionParser parser(option);
96                         map<string,string> parameters = parser.getParameters();
97                         
98                         ValidParameters validParameter;
99                         map<string,string>::iterator it;
100                         
101                         //check to make sure all parameters are valid for command
102                         for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) { 
103                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
104                         }
105                         
106                         //initialize outputTypes
107                         vector<string> tempOutNames;
108                         outputTypes["accnos"] = tempOutNames;
109                         
110                         //if the user changes the output directory command factory will send this info to us in the output parameter 
111                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = "";         }
112                         
113                         //if the user changes the input directory command factory will send this info to us in the output parameter 
114                         string inputDir = validParameter.validFile(parameters, "inputdir", false);              
115                         if (inputDir == "not found"){   inputDir = "";          }
116                         else {
117                                 string path;
118                                 it = parameters.find("alignreport");
119                                 //user has given a template file
120                                 if(it != parameters.end()){ 
121                                         path = m->hasPath(it->second);
122                                         //if the user has not given a path then, add inputdir. else leave path alone.
123                                         if (path == "") {       parameters["alignreport"] = inputDir + it->second;              }
124                                 }
125                                 
126                                 it = parameters.find("fasta");
127                                 //user has given a template file
128                                 if(it != parameters.end()){ 
129                                         path = m->hasPath(it->second);
130                                         //if the user has not given a path then, add inputdir. else leave path alone.
131                                         if (path == "") {       parameters["fasta"] = inputDir + it->second;            }
132                                 }
133                                 
134                                 it = parameters.find("list");
135                                 //user has given a template file
136                                 if(it != parameters.end()){ 
137                                         path = m->hasPath(it->second);
138                                         //if the user has not given a path then, add inputdir. else leave path alone.
139                                         if (path == "") {       parameters["list"] = inputDir + it->second;             }
140                                 }
141                                 
142                                 it = parameters.find("name");
143                                 //user has given a template file
144                                 if(it != parameters.end()){ 
145                                         path = m->hasPath(it->second);
146                                         //if the user has not given a path then, add inputdir. else leave path alone.
147                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
148                                 }
149                                 
150                                 it = parameters.find("group");
151                                 //user has given a template file
152                                 if(it != parameters.end()){ 
153                                         path = m->hasPath(it->second);
154                                         //if the user has not given a path then, add inputdir. else leave path alone.
155                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
156                                 }
157                                 
158                                 it = parameters.find("taxonomy");
159                                 //user has given a template file
160                                 if(it != parameters.end()){ 
161                                         path = m->hasPath(it->second);
162                                         //if the user has not given a path then, add inputdir. else leave path alone.
163                                         if (path == "") {       parameters["taxonomy"] = inputDir + it->second;         }
164                                 }
165                 
166                 it = parameters.find("count");
167                                 //user has given a template file
168                                 if(it != parameters.end()){ 
169                                         path = m->hasPath(it->second);
170                                         //if the user has not given a path then, add inputdir. else leave path alone.
171                                         if (path == "") {       parameters["count"] = inputDir + it->second;            }
172                                 }
173                 
174                 it = parameters.find("fastq");
175                                 //user has given a template file
176                                 if(it != parameters.end()){
177                                         path = m->hasPath(it->second);
178                                         //if the user has not given a path then, add inputdir. else leave path alone.
179                                         if (path == "") {       parameters["fastq"] = inputDir + it->second;            }
180                                 }
181                         }
182
183                         //check for required parameters
184                         fastafile = validParameter.validFile(parameters, "fasta", true);
185                         if (fastafile == "not open") { abort = true; }
186                         else if (fastafile == "not found") {  fastafile = "";  }
187                         else { m->setFastaFile(fastafile); }
188                         
189                         namefile = validParameter.validFile(parameters, "name", true);
190                         if (namefile == "not open") { abort = true; }
191                         else if (namefile == "not found") {  namefile = "";  }  
192                         else { m->setNameFile(namefile); }
193                         
194                         groupfile = validParameter.validFile(parameters, "group", true);
195                         if (groupfile == "not open") { abort = true; }
196                         else if (groupfile == "not found") {  groupfile = "";  }        
197                         else { m->setGroupFile(groupfile); }
198                         
199                         alignfile = validParameter.validFile(parameters, "alignreport", true);
200                         if (alignfile == "not open") { abort = true; }
201                         else if (alignfile == "not found") {  alignfile = "";  }
202                         
203                         listfile = validParameter.validFile(parameters, "list", true);
204                         if (listfile == "not open") { abort = true; }
205                         else if (listfile == "not found") {  listfile = "";  }
206                         else { m->setListFile(listfile); }
207                         
208                         taxfile = validParameter.validFile(parameters, "taxonomy", true);
209                         if (taxfile == "not open") { abort = true; }
210                         else if (taxfile == "not found") {  taxfile = "";  }
211                         else { m->setTaxonomyFile(taxfile); }
212             
213             countfile = validParameter.validFile(parameters, "count", true);
214                         if (countfile == "not open") { abort = true; }
215                         else if (countfile == "not found") {  countfile = "";  }
216                         else { m->setCountTableFile(countfile); }
217             
218             fastqfile = validParameter.validFile(parameters, "fastq", true);
219                         if (fastqfile == "not open") { abort = true; }
220                         else if (fastqfile == "not found") {  fastqfile = "";  }
221                         
222                         if ((fastqfile == "") && (countfile == "") && (fastafile == "") && (namefile == "") && (listfile == "") && (groupfile == "") && (alignfile == "") && (taxfile == ""))  { m->mothurOut("You must provide a file."); m->mothurOutEndLine(); abort = true; }
223                         
224                         int okay = 1;
225                         if (outputDir != "") { okay++; }
226                         if (inputDir != "") { okay++; }
227                         
228                         if (parameters.size() > okay) { m->mothurOut("You may only enter one file."); m->mothurOutEndLine(); abort = true;  }
229                 }
230
231         }
232         catch(exception& e) {
233                 m->errorOut(e, "ListSeqsCommand", "ListSeqsCommand");
234                 exit(1);
235         }
236 }
237 //**********************************************************************************************************************
238
239 int ListSeqsCommand::execute(){
240         try {
241                 
242                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
243                 
244                 //read functions fill names vector
245                 if (fastafile != "")            {       inputFileName = fastafile;      readFasta();    }
246         else if (fastqfile != "")       {       inputFileName = fastqfile;      readFastq();    }
247                 else if (namefile != "")        {       inputFileName = namefile;       readName();             }
248                 else if (groupfile != "")       {       inputFileName = groupfile;      readGroup();    }
249                 else if (alignfile != "")       {       inputFileName = alignfile;      readAlign();    }
250                 else if (listfile != "")        {       inputFileName = listfile;       readList();             }
251                 else if (taxfile != "")         {       inputFileName = taxfile;        readTax();              }
252         else if (countfile != "")       {       inputFileName = countfile;      readCount();    }
253                 
254                 if (m->control_pressed) { outputTypes.clear();  return 0; }
255                 
256                 //sort in alphabetical order
257                 sort(names.begin(), names.end());
258                 
259                 if (outputDir == "") {  outputDir += m->hasPath(inputFileName);  }
260                 
261         map<string, string> variables; 
262         variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(inputFileName));
263                 string outputFileName = getOutputFileName("accnos", variables);
264
265                 ofstream out;
266                 m->openOutputFile(outputFileName, out);
267                 outputNames.push_back(outputFileName); outputTypes["accnos"].push_back(outputFileName);
268                 
269                 //output to .accnos file
270                 for (int i = 0; i < names.size(); i++) {
271                         
272                         if (m->control_pressed) { outputTypes.clear(); out.close(); m->mothurRemove(outputFileName); return 0; }
273                         
274                         out << names[i] << endl;
275                 }
276                 out.close();
277                 
278                 if (m->control_pressed) { outputTypes.clear();  m->mothurRemove(outputFileName); return 0; }
279                 
280                 m->setAccnosFile(outputFileName);
281                 
282                 m->mothurOutEndLine();
283                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
284                 m->mothurOut(outputFileName); m->mothurOutEndLine();    
285                 m->mothurOutEndLine();
286                 
287                 //set accnos file as new current accnosfile
288                 string current = "";
289                 itTypes = outputTypes.find("accnos");
290                 if (itTypes != outputTypes.end()) {
291                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setAccnosFile(current); }
292                 }
293                 
294                 return 0;               
295         }
296
297         catch(exception& e) {
298                 m->errorOut(e, "ListSeqsCommand", "execute");
299                 exit(1);
300         }
301 }
302 //**********************************************************************************************************************
303 int ListSeqsCommand::readFastq(){
304         try {
305                 
306                 ifstream in;
307                 m->openInputFile(fastqfile, in);
308                 string name;
309                 
310                 //ofstream out;
311                 //string newFastaName = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "numsAdded.fasta";
312                 //m->openOutputFile(newFastaName, out);
313                 int count = 1;
314                 //string lastName = "";
315                 
316                 while(!in.eof()){
317                         
318                         if (m->control_pressed) { in.close(); return 0; }
319                         
320                         //read sequence name
321                         string name = m->getline(in); m->gobble(in);
322                         
323                         if (name[0] == '@') {
324                 vector<string> splits = m->splitWhiteSpace(name);
325                 name = splits[0];
326                 name = name.substr(1);
327                 m->checkName(name);
328                 names.push_back(name);
329                 //get rest of lines
330                 name = m->getline(in); m->gobble(in);
331                 name = m->getline(in); m->gobble(in);
332                 name = m->getline(in); m->gobble(in);
333             }
334                         
335                         m->gobble(in);
336                         if (m->debug) { count++; cout << "[DEBUG]: count = " + toString(count) + ", name = " + name + "\n"; }
337                 }
338                 in.close();
339                 //out.close();
340                 
341                 return 0;
342         
343         }
344         catch(exception& e) {
345                 m->errorOut(e, "ListSeqsCommand", "readFastq");
346                 exit(1);
347         }
348 }
349
350 //**********************************************************************************************************************
351 int ListSeqsCommand::readFasta(){
352         try {
353                 
354                 ifstream in;
355                 m->openInputFile(fastafile, in);
356                 string name;
357                 
358                 //ofstream out;
359                 //string newFastaName = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "numsAdded.fasta";
360                 //m->openOutputFile(newFastaName, out);
361                 int count = 1;
362                 //string lastName = "";
363                 
364                 while(!in.eof()){
365                         
366                         if (m->control_pressed) { in.close(); return 0; }
367                         
368                         Sequence currSeq(in);
369                         name = currSeq.getName();
370                         
371                         if (name != "") {  names.push_back(name);  }
372                         
373                         m->gobble(in);
374                         if (m->debug) { count++; cout << "[DEBUG]: count = " + toString(count) + ", name = " + currSeq.getName() + "\n"; }
375                 }
376                 in.close();     
377                 //out.close();
378                 
379                 return 0;
380
381         }
382         catch(exception& e) {
383                 m->errorOut(e, "ListSeqsCommand", "readFasta");
384                 exit(1);
385         }
386 }
387 //**********************************************************************************************************************
388 int ListSeqsCommand::readList(){
389         try {
390                 ifstream in;
391                 m->openInputFile(listfile, in);
392                 
393                 if(!in.eof()){
394                         //read in list vector
395                         ListVector list(in);
396                         
397                         //for each bin
398                         for (int i = 0; i < list.getNumBins(); i++) {
399                                 string binnames = list.get(i);
400                                 
401                                 if (m->control_pressed) { in.close(); return 0; }
402                                 
403                                 m->splitAtComma(binnames, names);
404                         }
405                 }
406                 in.close();     
407                 
408                 return 0;
409                 
410         }
411         catch(exception& e) {
412                 m->errorOut(e, "ListSeqsCommand", "readList");
413                 exit(1);
414         }
415 }
416
417 //**********************************************************************************************************************
418 int ListSeqsCommand::readName(){
419         try {
420                 
421                 ifstream in;
422                 m->openInputFile(namefile, in);
423                 string name, firstCol, secondCol;
424                 
425                 while(!in.eof()){
426                 
427                         if (m->control_pressed) { in.close(); return 0; }
428
429                         in >> firstCol;                         
430                         in >> secondCol;                        
431                         
432                         //parse second column saving each name
433                         m->splitAtComma(secondCol, names);
434                         
435                         m->gobble(in);
436                 }
437                 in.close();
438                 return 0;
439                 
440         }
441         catch(exception& e) {
442                 m->errorOut(e, "ListSeqsCommand", "readName");
443                 exit(1);
444         }
445 }
446
447 //**********************************************************************************************************************
448 int ListSeqsCommand::readGroup(){
449         try {
450         
451                 ifstream in;
452                 m->openInputFile(groupfile, in);
453                 string name, group;
454                 
455                 while(!in.eof()){
456                         
457                         if (m->control_pressed) { in.close(); return 0; }
458                         
459                         in >> name;     m->gobble(in);                  //read from first column
460                         in >> group;                    //read from second column
461                         
462                         names.push_back(name);
463                                         
464                         m->gobble(in);
465                 }
466                 in.close();
467                 return 0;
468
469         }
470         catch(exception& e) {
471                 m->errorOut(e, "ListSeqsCommand", "readGroup");
472                 exit(1);
473         }
474 }
475 //**********************************************************************************************************************
476 int ListSeqsCommand::readCount(){
477         try {
478                 CountTable ct;
479                 ct.readTable(countfile, false, false);
480         
481         if (m->control_pressed) { return 0; }
482         
483         names = ct.getNamesOfSeqs();
484         
485         return 0;
486         
487         }
488         catch(exception& e) {
489                 m->errorOut(e, "ListSeqsCommand", "readCount");
490                 exit(1);
491         }
492 }
493 //**********************************************************************************************************************
494 //alignreport file has a column header line then all other lines contain 16 columns.  we just want the first column since that contains the name
495 int ListSeqsCommand::readAlign(){
496         try {
497         
498                 ifstream in;
499                 m->openInputFile(alignfile, in);
500                 string name, junk;
501                 
502                 //read column headers
503                 for (int i = 0; i < 16; i++) {  
504                         if (!in.eof())  {       in >> junk;             }
505                         else                    {       break;                  }
506                 }
507                 //m->getline(in);
508                 
509                 while(!in.eof()){
510                 
511                         if (m->control_pressed) { in.close(); return 0; }
512
513                         in >> name;                             //read from first column
514                         //m->getline(in);
515                         //read rest
516                         for (int i = 0; i < 15; i++) {  
517                                 if (!in.eof())  {       in >> junk;             }
518                                 else                    {       break;                  }
519                         }
520                         
521                         names.push_back(name);
522                                         
523                         m->gobble(in);
524                 }
525                 in.close();
526                 
527                 return 0;
528
529                 
530         }
531         catch(exception& e) {
532                 m->errorOut(e, "ListSeqsCommand", "readAlign");
533                 exit(1);
534         }
535 }
536 //**********************************************************************************************************************
537 int ListSeqsCommand::readTax(){
538         try {
539                 
540                 ifstream in;
541                 m->openInputFile(taxfile, in);
542                 string name, firstCol, secondCol;
543                 
544                 while(!in.eof()){
545                 
546                         if (m->control_pressed) { in.close(); return 0; }
547
548                         in >> firstCol;                         
549                         in >> secondCol;                        
550                         
551                         names.push_back(firstCol);
552                         
553                         m->gobble(in);
554                         
555                 }
556                 in.close();
557                 
558                 return 0;
559                 
560         }
561         catch(exception& e) {
562                 m->errorOut(e, "ListSeqsCommand", "readTax");
563                 exit(1);
564         }
565 }
566 //**********************************************************************************************************************