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 string CountSeqsCommand::getOutputFileNameTag(string type, string inputName=""){
  51         try {
  52         string outputFileName = "";
  53                 map<string, vector<string> >::iterator it;
  54
  55         //is this a type this command creates
  56         it = outputTypes.find(type);
  57         if (it == outputTypes.end()) {  m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
  58         else {
  59             if (type == "summary") {  outputFileName =  "seq.count"; }
  60             else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true;  }
  61         }
  62         return outputFileName;
  63         }
  64         catch(exception& e) {
  65                 m->errorOut(e, "CountSeqsCommand", "getOutputFileNameTag");
  66                 exit(1);
  67         }
  68 }
  69 //**********************************************************************************************************************
  70 CountSeqsCommand::CountSeqsCommand(){
  71         try {
  72                 abort = true; calledHelp = true;
  73                 setParameters();
  74                 vector<string> tempOutNames;
  75                 outputTypes["summary"] = tempOutNames;
  76         }
  77         catch(exception& e) {
  78                 m->errorOut(e, "CountSeqsCommand", "CountSeqsCommand");
  79                 exit(1);
  80         }
  81 }
  82 //**********************************************************************************************************************
  83
  84 CountSeqsCommand::CountSeqsCommand(string option)  {
  85         try {
  86                 abort = false; calledHelp = false;
  87
  88                 //allow user to run help
  89                 if(option == "help") { help(); abort = true; calledHelp = true; }
  90                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  91                 else {
  92                         vector<string> myArray = setParameters();
  93
  94                         OptionParser parser(option);
  95                         map<string,string> parameters = parser.getParameters();
  96
  97                         ValidParameters validParameter;
  98                         map<string,string>::iterator it;
  99
 100                         //check to make sure all parameters are valid for command
 101                         for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
 102                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 103                         }
 104
 105                         //initialize outputTypes
 106                         vector<string> tempOutNames;
 107                         outputTypes["summary"] = tempOutNames;
 108
 109
 110                         //if the user changes the input directory command factory will send this info to us in the output parameter
 111                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 112                         if (inputDir == "not found"){   inputDir = "";          }
 113                         else {
 114                                 string path;
 115                                 it = parameters.find("name");
 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["name"] = inputDir + it->second;             }
 121                                 }
 122
 123                                 it = parameters.find("group");
 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["group"] = inputDir + it->second;            }
 129                                 }
 130                         }
 131
 132                         //check for required parameters
 133                         namefile = validParameter.validFile(parameters, "name", true);
 134                         if (namefile == "not open") { namefile = ""; abort = true; }
 135                         else if (namefile == "not found"){
 136                                 namefile = m->getNameFile();
 137                                 if (namefile != "") { m->mothurOut("Using " + namefile + " as input file for the name parameter."); m->mothurOutEndLine(); }
 138                                 else {  m->mothurOut("You have no current namefile and the name parameter is required."); m->mothurOutEndLine(); abort = true; }
 139                         }else { m->setNameFile(namefile); }
 140
 141                         groupfile = validParameter.validFile(parameters, "group", true);
 142                         if (groupfile == "not open") { abort = true; }
 143                         else if (groupfile == "not found") {  groupfile = "";  }
 144                         else { m->setGroupFile(groupfile); }
 145
 146                         groups = validParameter.validFile(parameters, "groups", false);
 147                         if (groups == "not found") { groups = "all"; }
 148                         m->splitAtDash(groups, Groups);
 149
 150                         //if the user changes the output directory command factory will send this info to us in the output parameter
 151                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(namefile);               }
 152
 153                 }
 154
 155         }
 156         catch(exception& e) {
 157                 m->errorOut(e, "CountSeqsCommand", "CountSeqsCommand");
 158                 exit(1);
 159         }
 160 }
 161 //**********************************************************************************************************************
 162
 163 int CountSeqsCommand::execute(){
 164         try {
 165
 166                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 167
 168                 ofstream out;
 169                 string outputFileName = outputDir + m->getRootName(m->getSimpleName(namefile)) + getOutputFileNameTag("summary");
 170                 m->openOutputFile(outputFileName, out); outputTypes["summary"].push_back(outputFileName);
 171                 out << "Representative_Sequence\ttotal\t";
 172
 173                 GroupMap* groupMap;
 174                 if (groupfile != "") {
 175                         groupMap = new GroupMap(groupfile); groupMap->readMap();
 176
 177                         //make sure groups are valid. takes care of user setting groupNames that are invalid or setting groups=all
 178                         SharedUtil* util = new SharedUtil();
 179                         vector<string> nameGroups = groupMap->getNamesOfGroups();
 180                         util->setGroups(Groups, nameGroups);
 181                         delete util;
 182
 183                         //sort groupNames so that the group title match the counts below, this is needed because the map object automatically sorts
 184                         sort(Groups.begin(), Groups.end());
 185
 186                         //print groupNames
 187                         for (int i = 0; i < Groups.size(); i++) {
 188                                 out << Groups[i] << '\t';
 189                         }
 190                 }
 191                 out << endl;
 192
 193                 //open input file
 194                 ifstream in;
 195                 m->openInputFile(namefile, in);
 196
 197                 int total = 0;
 198                 while (!in.eof()) {
 199                         if (m->control_pressed) { break; }
 200
 201                         string firstCol, secondCol;
 202                         in >> firstCol >> secondCol; m->gobble(in);
 203
 204                         vector<string> names;
 205                         m->splitAtChar(secondCol, names, ',');
 206
 207                         if (groupfile != "") {
 208                                 //set to 0
 209                                 map<string, int> groupCounts;
 210                                 int total = 0;
 211                                 for (int i = 0; i < Groups.size(); i++) { groupCounts[Groups[i]] = 0; }
 212
 213                                 //get counts for each of the users groups
 214                                 for (int i = 0; i < names.size(); i++) {
 215                                         string group = groupMap->getGroup(names[i]);
 216
 217                                         if (group == "not found") { m->mothurOut("[ERROR]: " + names[i] + " is not in your groupfile, please correct."); m->mothurOutEndLine(); }
 218                                         else {
 219                                                 map<string, int>::iterator it = groupCounts.find(group);
 220
 221                                                 //if not found, then this sequence is not from a group we care about
 222                                                 if (it != groupCounts.end()) {
 223                                                         it->second++;
 224                                                         total++;
 225                                                 }
 226                                         }
 227                                 }
 228
 229                                 if (total != 0) {
 230                                         out << firstCol << '\t' << total << '\t';
 231                                         for (map<string, int>::iterator it = groupCounts.begin(); it != groupCounts.end(); it++) {
 232                                                 out << it->second << '\t';
 233                                         }
 234                                         out << endl;
 235                                 }
 236                         }else {
 237                                 out << firstCol << '\t' << names.size() << endl;
 238                         }
 239
 240                         total += names.size();
 241                 }
 242                 in.close();
 243
 244                 if (groupfile != "") { delete groupMap; }
 245
 246                 if (m->control_pressed) { m->mothurRemove(outputFileName); return 0; }
 247
 248         m->mothurOutEndLine();
 249                 m->mothurOut("Total number of sequences: " + toString(total)); m->mothurOutEndLine();
 250                 m->mothurOutEndLine();
 251                 m->mothurOut("Output File Name: "); m->mothurOutEndLine();
 252                 m->mothurOut(outputFileName); m->mothurOutEndLine();
 253                 m->mothurOutEndLine();
 254
 255                 return 0;
 256         }
 257
 258         catch(exception& e) {
 259                 m->errorOut(e, "CountSeqsCommand", "execute");
 260                 exit(1);
 261         }
 262 }
 263 //**********************************************************************************************************************