deconvolutecommand.cpp

   1 /*
   2  *  deconvolute.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 1/21/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "deconvolutecommand.h"
  11 #include "sequence.hpp"
  12
  13 //**********************************************************************************************************************
  14 vector<string> DeconvoluteCommand::setParameters(){
  15         try {
  16                 CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pfasta);
  17                 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
  18                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  19                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  20
  21                 vector<string> myArray;
  22                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  23                 return myArray;
  24         }
  25         catch(exception& e) {
  26                 m->errorOut(e, "DeconvoluteCommand", "setParameters");
  27                 exit(1);
  28         }
  29 }
  30 //**********************************************************************************************************************
  31 string DeconvoluteCommand::getHelpString(){
  32         try {
  33                 string helpString = "";
  34                 helpString += "The unique.seqs command reads a fastafile and creates a namesfile.\n";
  35                 helpString += "It creates a file where the first column is the groupname and the second column is a list of sequence names who have the same sequence. \n";
  36                 helpString += "If the sequence is unique the second column will just contain its name. \n";
  37                 helpString += "The unique.seqs command parameters are fasta and name.  fasta is required, unless there is a valid current fasta file.\n";
  38                 helpString += "The unique.seqs command should be in the following format: \n";
  39                 helpString += "unique.seqs(fasta=yourFastaFile) \n";
  40                 return helpString;
  41         }
  42         catch(exception& e) {
  43                 m->errorOut(e, "DeconvoluteCommand", "getHelpString");
  44                 exit(1);
  45         }
  46 }
  47 //**********************************************************************************************************************
  48 string DeconvoluteCommand::getOutputFileNameTag(string type, string inputName=""){
  49         try {
  50         string outputFileName = "";
  51                 map<string, vector<string> >::iterator it;
  52
  53         //is this a type this command creates
  54         it = outputTypes.find(type);
  55         if (it == outputTypes.end()) {  m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
  56         else {
  57             if (type == "fasta") {  outputFileName =  "unique" + m->getExtension(inputName); }
  58             else if (type == "name") {  outputFileName =  "names"; }
  59             else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true;  }
  60         }
  61         return outputFileName;
  62         }
  63         catch(exception& e) {
  64                 m->errorOut(e, "DeconvoluteCommand", "getOutputFileNameTag");
  65                 exit(1);
  66         }
  67 }
  68 //**********************************************************************************************************************
  69 DeconvoluteCommand::DeconvoluteCommand(){
  70         try {
  71                 abort = true; calledHelp = true;
  72                 setParameters();
  73                 vector<string> tempOutNames;
  74                 outputTypes["fasta"] = tempOutNames;
  75                 outputTypes["name"] = tempOutNames;
  76         }
  77         catch(exception& e) {
  78                 m->errorOut(e, "DeconvoluteCommand", "DeconvoluteCommand");
  79                 exit(1);
  80         }
  81 }
  82 /**************************************************************************************/
  83 DeconvoluteCommand::DeconvoluteCommand(string option)  {
  84         try {
  85                 abort = false; calledHelp = false;
  86
  87                 //allow user to run help
  88                 if(option == "help") { help(); abort = true; calledHelp = true; }
  89                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  90
  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 (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["fasta"] = tempOutNames;
 108                         outputTypes["name"] = tempOutNames;
 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("fasta");
 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["fasta"] = inputDir + it->second;            }
 121                                 }
 122
 123                                 it = parameters.find("name");
 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["name"] = inputDir + it->second;             }
 129                                 }
 130                         }
 131
 132
 133                         //check for required parameters
 134                         inFastaName = validParameter.validFile(parameters, "fasta", true);
 135                         if (inFastaName == "not open") { abort = true; }
 136                         else if (inFastaName == "not found") {
 137                                 inFastaName = m->getFastaFile();
 138                                 if (inFastaName != "") { m->mothurOut("Using " + inFastaName + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
 139                                 else {  m->mothurOut("You have no current fastafile and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
 140                         }else { m->setFastaFile(inFastaName); }
 141
 142                         //if the user changes the output directory command factory will send this info to us in the output parameter
 143                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
 144                                 outputDir = "";
 145                                 outputDir += m->hasPath(inFastaName); //if user entered a file with a path then preserve it
 146                         }
 147
 148                         oldNameMapFName = validParameter.validFile(parameters, "name", true);
 149                         if (oldNameMapFName == "not open") { oldNameMapFName = ""; abort = true; }
 150                         else if (oldNameMapFName == "not found"){       oldNameMapFName = "";   }
 151                         else { m->setNameFile(oldNameMapFName); }
 152
 153                         if (oldNameMapFName == "") {
 154                                 vector<string> files; files.push_back(inFastaName);
 155                                 parser.getNameFile(files);
 156                         }
 157
 158                 }
 159
 160         }
 161         catch(exception& e) {
 162                 m->errorOut(e, "DeconvoluteCommand", "DeconvoluteCommand");
 163                 exit(1);
 164         }
 165 }
 166 /**************************************************************************************/
 167 int DeconvoluteCommand::execute() {
 168         try {
 169
 170                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 171
 172                 //prepare filenames and open files
 173                 string outNameFile = outputDir + m->getRootName(m->getSimpleName(inFastaName)) + getOutputFileNameTag("name");
 174                 string outFastaFile = outputDir + m->getRootName(m->getSimpleName(inFastaName)) + getOutputFileNameTag("fasta", inFastaName);
 175
 176                 map<string, string> nameMap;
 177                 map<string, string>::iterator itNames;
 178                 if (oldNameMapFName != "")  {
 179             m->readNames(oldNameMapFName, nameMap);
 180             if (oldNameMapFName == outNameFile){ outNameFile = outputDir + m->getRootName(m->getSimpleName(inFastaName)) + "unique." + getOutputFileNameTag("name");   }
 181         }
 182
 183                 if (m->control_pressed) { return 0; }
 184
 185                 ifstream in;
 186                 m->openInputFile(inFastaName, in);
 187
 188                 ofstream outFasta;
 189                 m->openOutputFile(outFastaFile, outFasta);
 190
 191                 map<string, string> sequenceStrings; //sequenceString -> list of names.  "atgc...." -> seq1,seq2,seq3.
 192                 map<string, string>::iterator itStrings;
 193                 set<string> nameInFastaFile; //for sanity checking
 194                 set<string>::iterator itname;
 195                 vector<string> nameFileOrder;
 196                 int count = 0;
 197                 while (!in.eof()) {
 198
 199                         if (m->control_pressed) { in.close(); outFasta.close(); m->mothurRemove(outFastaFile); return 0; }
 200
 201                         Sequence seq(in);
 202
 203                         if (seq.getName() != "") {
 204
 205                                 //sanity checks
 206                                 itname = nameInFastaFile.find(seq.getName());
 207                                 if (itname == nameInFastaFile.end()) { nameInFastaFile.insert(seq.getName());  }
 208                                 else { m->mothurOut("[ERROR]: You already have a sequence named " + seq.getName() + " in your fasta file, sequence names must be unique, please correct."); m->mothurOutEndLine(); }
 209
 210                                 itStrings = sequenceStrings.find(seq.getAligned());
 211
 212                                 if (itStrings == sequenceStrings.end()) { //this is a new unique sequence
 213                                         //output to unique fasta file
 214                                         seq.printSequence(outFasta);
 215
 216                                         if (oldNameMapFName != "") {
 217                                                 itNames = nameMap.find(seq.getName());
 218
 219                                                 if (itNames == nameMap.end()) { //namefile and fastafile do not match
 220                                                         m->mothurOut("[ERROR]: " + seq.getName() + " is in your fasta file, and not in your namefile, please correct."); m->mothurOutEndLine();
 221                                                 }else {
 222                                                         sequenceStrings[seq.getAligned()] = itNames->second;
 223                                                         nameFileOrder.push_back(seq.getAligned());
 224                                                 }
 225                                         }else { sequenceStrings[seq.getAligned()] = seq.getName();      nameFileOrder.push_back(seq.getAligned()); }
 226                                 }else { //this is a dup
 227                                         if (oldNameMapFName != "") {
 228                                                 itNames = nameMap.find(seq.getName());
 229
 230                                                 if (itNames == nameMap.end()) { //namefile and fastafile do not match
 231                                                         m->mothurOut("[ERROR]: " + seq.getName() + " is in your fasta file, and not in your namefile, please correct."); m->mothurOutEndLine();
 232                                                 }else {
 233                                                         sequenceStrings[seq.getAligned()] += "," + itNames->second;
 234                                                 }
 235                                         }else { sequenceStrings[seq.getAligned()] += "," + seq.getName();       }
 236                                 }
 237
 238                                 count++;
 239                         }
 240
 241                         m->gobble(in);
 242
 243                         if(count % 1000 == 0)   { m->mothurOut(toString(count) + "\t" + toString(sequenceStrings.size())); m->mothurOutEndLine();       }
 244                 }
 245
 246                 if(count % 1000 != 0)   { m->mothurOut(toString(count) + "\t" + toString(sequenceStrings.size())); m->mothurOutEndLine();       }
 247
 248                 in.close();
 249                 outFasta.close();
 250
 251                 if (m->control_pressed) { m->mothurRemove(outFastaFile); return 0; }
 252
 253                 //print new names file
 254                 ofstream outNames;
 255                 m->openOutputFile(outNameFile, outNames);
 256
 257                 for (int i = 0; i < nameFileOrder.size(); i++) {
 258                 //for (itStrings = sequenceStrings.begin(); itStrings != sequenceStrings.end(); itStrings++) {
 259                         if (m->control_pressed) { outputTypes.clear(); m->mothurRemove(outFastaFile); outNames.close(); m->mothurRemove(outNameFile); return 0; }
 260
 261                         itStrings = sequenceStrings.find(nameFileOrder[i]);
 262
 263                         if (itStrings != sequenceStrings.end()) {
 264                                 //get rep name
 265                                 int pos = (itStrings->second).find_first_of(',');
 266
 267                                 if (pos == string::npos) { // only reps itself
 268                                         outNames << itStrings->second << '\t' << itStrings->second << endl;
 269                                 }else {
 270                                         outNames << (itStrings->second).substr(0, pos) << '\t' << itStrings->second << endl;
 271                                 }
 272                         }else{ m->mothurOut("[ERROR]: mismatch in namefile print."); m->mothurOutEndLine(); m->control_pressed = true; }
 273                 }
 274                 outNames.close();
 275
 276                 if (m->control_pressed) { outputTypes.clear(); m->mothurRemove(outFastaFile); m->mothurRemove(outNameFile); return 0; }
 277
 278                 m->mothurOutEndLine();
 279                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 280                 m->mothurOut(outFastaFile); m->mothurOutEndLine();
 281                 m->mothurOut(outNameFile); m->mothurOutEndLine();
 282                 outputNames.push_back(outFastaFile);  outputNames.push_back(outNameFile); outputTypes["fasta"].push_back(outFastaFile);  outputTypes["name"].push_back(outNameFile);
 283                 m->mothurOutEndLine();
 284
 285                 //set fasta file as new current fastafile
 286                 string current = "";
 287                 itTypes = outputTypes.find("fasta");
 288                 if (itTypes != outputTypes.end()) {
 289                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
 290                 }
 291
 292                 itTypes = outputTypes.find("name");
 293                 if (itTypes != outputTypes.end()) {
 294                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setNameFile(current); }
 295                 }
 296
 297                 return 0;
 298         }
 299         catch(exception& e) {
 300                 m->errorOut(e, "DeconvoluteCommand", "execute");
 301                 exit(1);
 302         }
 303 }
 304 /**************************************************************************************/