countseqscommand.cpp

   1 /*
   2  *  countseqscommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 6/1/11.
   6  *  Copyright 2011 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "countseqscommand.h"
  11 #include "groupmap.h"
  12 #include "sharedutilities.h"
  13
  14 //**********************************************************************************************************************
  15 vector<string> CountSeqsCommand::setParameters(){
  16         try {
  17                 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pname);
  18                 CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pgroup);
  19                 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
  20                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  21                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  22
  23                 vector<string> myArray;
  24                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  25                 return myArray;
  26         }
  27         catch(exception& e) {
  28                 m->errorOut(e, "CountSeqsCommand", "setParameters");
  29                 exit(1);
  30         }
  31 }
  32 //**********************************************************************************************************************
  33 string CountSeqsCommand::getHelpString(){
  34         try {
  35                 string helpString = "";
  36                 helpString += "The count.seqs command reads a name file and outputs a .seq.count file.  You may also provide a group file to get the counts broken down by group.\n";
  37                 helpString += "The groups parameter allows you to indicate which groups you want to include in the counts, by default all groups in your groupfile are used.\n";
  38                 helpString += "When you use the groups parameter and a sequence does not represent any sequences from the groups you specify it is not included in the .count.summary file.\n";
  39                 helpString += "The count.seqs command should be in the following format: count.seqs(name=yourNameFile).\n";
  40                 helpString += "Example count.seqs(name=amazon.names).\n";
  41                 helpString += "Note: No spaces between parameter labels (i.e. name), '=' and parameters (i.e.yourNameFile).\n";
  42                 return helpString;
  43         }
  44         catch(exception& e) {
  45                 m->errorOut(e, "CountSeqsCommand", "getHelpString");
  46                 exit(1);
  47         }
  48 }
  49
  50 //**********************************************************************************************************************
  51 CountSeqsCommand::CountSeqsCommand(){
  52         try {
  53                 abort = true; calledHelp = true;
  54                 setParameters();
  55                 vector<string> tempOutNames;
  56                 outputTypes["summary"] = tempOutNames;
  57         }
  58         catch(exception& e) {
  59                 m->errorOut(e, "CountSeqsCommand", "CountSeqsCommand");
  60                 exit(1);
  61         }
  62 }
  63 //**********************************************************************************************************************
  64
  65 CountSeqsCommand::CountSeqsCommand(string option)  {
  66         try {
  67                 abort = false; calledHelp = false;
  68
  69                 //allow user to run help
  70                 if(option == "help") { help(); abort = true; calledHelp = true; }
  71                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  72                 else {
  73                         vector<string> myArray = setParameters();
  74
  75                         OptionParser parser(option);
  76                         map<string,string> parameters = parser.getParameters();
  77
  78                         ValidParameters validParameter;
  79                         map<string,string>::iterator it;
  80
  81                         //check to make sure all parameters are valid for command
  82                         for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
  83                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  84                         }
  85
  86                         //initialize outputTypes
  87                         vector<string> tempOutNames;
  88                         outputTypes["summary"] = tempOutNames;
  89
  90
  91                         //if the user changes the input directory command factory will send this info to us in the output parameter
  92                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  93                         if (inputDir == "not found"){   inputDir = "";          }
  94                         else {
  95                                 string path;
  96                                 it = parameters.find("name");
  97                                 //user has given a template file
  98                                 if(it != parameters.end()){
  99                                         path = m->hasPath(it->second);
 100                                         //if the user has not given a path then, add inputdir. else leave path alone.
 101                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
 102                                 }
 103
 104                                 it = parameters.find("group");
 105                                 //user has given a template file
 106                                 if(it != parameters.end()){
 107                                         path = m->hasPath(it->second);
 108                                         //if the user has not given a path then, add inputdir. else leave path alone.
 109                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
 110                                 }
 111                         }
 112
 113                         //check for required parameters
 114                         namefile = validParameter.validFile(parameters, "name", true);
 115                         if (namefile == "not open") { namefile = ""; abort = true; }
 116                         else if (namefile == "not found"){
 117                                 namefile = m->getNameFile();
 118                                 if (namefile != "") { m->mothurOut("Using " + namefile + " as input file for the name parameter."); m->mothurOutEndLine(); }
 119                                 else {  m->mothurOut("You have no current namefile and the name parameter is required."); m->mothurOutEndLine(); abort = true; }
 120                         }else { m->setNameFile(namefile); }
 121
 122                         groupfile = validParameter.validFile(parameters, "group", true);
 123                         if (groupfile == "not open") { abort = true; }
 124                         else if (groupfile == "not found") {  groupfile = "";  }
 125                         else { m->setGroupFile(groupfile); }
 126
 127                         groups = validParameter.validFile(parameters, "groups", false);
 128                         if (groups == "not found") { groups = "all"; }
 129                         m->splitAtDash(groups, Groups);
 130
 131                         //if the user changes the output directory command factory will send this info to us in the output parameter
 132                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(namefile);               }
 133
 134                 }
 135
 136         }
 137         catch(exception& e) {
 138                 m->errorOut(e, "CountSeqsCommand", "CountSeqsCommand");
 139                 exit(1);
 140         }
 141 }
 142 //**********************************************************************************************************************
 143
 144 int CountSeqsCommand::execute(){
 145         try {
 146
 147                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 148
 149                 ofstream out;
 150                 string outputFileName = outputDir + m->getRootName(m->getSimpleName(namefile)) + "seq.count";
 151                 m->openOutputFile(outputFileName, out); outputTypes["summary"].push_back(outputFileName);
 152                 out << "Representative_Sequence\ttotal\t";
 153
 154                 GroupMap* groupMap;
 155                 if (groupfile != "") {
 156                         groupMap = new GroupMap(groupfile); groupMap->readMap();
 157
 158                         //make sure groups are valid. takes care of user setting groupNames that are invalid or setting groups=all
 159                         SharedUtil* util = new SharedUtil();
 160                         vector<string> nameGroups = groupMap->getNamesOfGroups();
 161                         util->setGroups(Groups, nameGroups);
 162                         delete util;
 163
 164                         //sort groupNames so that the group title match the counts below, this is needed because the map object automatically sorts
 165                         sort(Groups.begin(), Groups.end());
 166
 167                         //print groupNames
 168                         for (int i = 0; i < Groups.size(); i++) {
 169                                 out << Groups[i] << '\t';
 170                         }
 171                 }
 172                 out << endl;
 173
 174                 //open input file
 175                 ifstream in;
 176                 m->openInputFile(namefile, in);
 177
 178                 int total = 0;
 179                 while (!in.eof()) {
 180                         if (m->control_pressed) { break; }
 181
 182                         string firstCol, secondCol;
 183                         in >> firstCol >> secondCol; m->gobble(in);
 184
 185                         vector<string> names;
 186                         m->splitAtChar(secondCol, names, ',');
 187
 188                         if (groupfile != "") {
 189                                 //set to 0
 190                                 map<string, int> groupCounts;
 191                                 int total = 0;
 192                                 for (int i = 0; i < Groups.size(); i++) { groupCounts[Groups[i]] = 0; }
 193
 194                                 //get counts for each of the users groups
 195                                 for (int i = 0; i < names.size(); i++) {
 196                                         string group = groupMap->getGroup(names[i]);
 197
 198                                         if (group == "not found") { m->mothurOut("[ERROR]: " + names[i] + " is not in your groupfile, please correct."); m->mothurOutEndLine(); }
 199                                         else {
 200                                                 map<string, int>::iterator it = groupCounts.find(group);
 201
 202                                                 //if not found, then this sequence is not from a group we care about
 203                                                 if (it != groupCounts.end()) {
 204                                                         it->second++;
 205                                                         total++;
 206                                                 }
 207                                         }
 208                                 }
 209
 210                                 if (total != 0) {
 211                                         out << firstCol << '\t' << total << '\t';
 212                                         for (map<string, int>::iterator it = groupCounts.begin(); it != groupCounts.end(); it++) {
 213                                                 out << it->second << '\t';
 214                                         }
 215                                         out << endl;
 216                                 }
 217                         }else {
 218                                 out << firstCol << '\t' << names.size() << endl;
 219                         }
 220
 221                         total += names.size();
 222                 }
 223                 in.close();
 224
 225                 if (groupfile != "") { delete groupMap; }
 226
 227                 if (m->control_pressed) { m->mothurRemove(outputFileName); return 0; }
 228
 229         m->mothurOutEndLine();
 230                 m->mothurOut("Total number of sequences: " + toString(total)); m->mothurOutEndLine();
 231                 m->mothurOutEndLine();
 232                 m->mothurOut("Output File Name: "); m->mothurOutEndLine();
 233                 m->mothurOut(outputFileName); m->mothurOutEndLine();
 234                 m->mothurOutEndLine();
 235
 236                 return 0;
 237         }
 238
 239         catch(exception& e) {
 240                 m->errorOut(e, "CountSeqsCommand", "execute");
 241                 exit(1);
 242         }
 243 }
 244 //**********************************************************************************************************************