chopseqscommand.cpp

   1 /*
   2  *  chopseqscommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 5/10/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "chopseqscommand.h"
  11 #include "sequence.hpp"
  12
  13 //**********************************************************************************************************************
  14 vector<string> ChopSeqsCommand::getValidParameters(){
  15         try {
  16                 string AlignArray[] =  {"fasta","short","numbases","countgaps","keep","outputdir","inputdir"};
  17                 vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
  18                 return myArray;
  19         }
  20         catch(exception& e) {
  21                 m->errorOut(e, "ChopSeqsCommand", "getValidParameters");
  22                 exit(1);
  23         }
  24 }
  25 //**********************************************************************************************************************
  26 ChopSeqsCommand::ChopSeqsCommand(){
  27         try {
  28                 abort = true; calledHelp = true;
  29                 vector<string> tempOutNames;
  30                 outputTypes["fasta"] = tempOutNames;
  31                 outputTypes["accnos"] = tempOutNames;
  32         }
  33         catch(exception& e) {
  34                 m->errorOut(e, "ChopSeqsCommand", "ChopSeqsCommand");
  35                 exit(1);
  36         }
  37 }
  38 //**********************************************************************************************************************
  39 vector<string> ChopSeqsCommand::getRequiredParameters(){
  40         try {
  41                 string Array[] =  {"fasta","numbases"};
  42                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  43                 return myArray;
  44         }
  45         catch(exception& e) {
  46                 m->errorOut(e, "ChopSeqsCommand", "getRequiredParameters");
  47                 exit(1);
  48         }
  49 }
  50 //**********************************************************************************************************************
  51 vector<string> ChopSeqsCommand::getRequiredFiles(){
  52         try {
  53                 vector<string> myArray;
  54                 return myArray;
  55         }
  56         catch(exception& e) {
  57                 m->errorOut(e, "ChopSeqsCommand", "getRequiredFiles");
  58                 exit(1);
  59         }
  60 }
  61 //**********************************************************************************************************************
  62 ChopSeqsCommand::ChopSeqsCommand(string option)  {
  63         try {
  64                 abort = false; calledHelp = false;
  65
  66                 //allow user to run help
  67                 if(option == "help") { help(); abort = true; calledHelp = true; }
  68
  69                 else {
  70                         //valid paramters for this command
  71                         string Array[] =  {"fasta","numbases","countgaps","keep","short","outputdir","inputdir"};
  72                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  73
  74                         OptionParser parser(option);
  75                         map<string,string> parameters = parser.getParameters();
  76
  77                         ValidParameters validParameter;
  78                         map<string,string>::iterator it;
  79
  80                         //check to make sure all parameters are valid for command
  81                         for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
  82                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  83                         }
  84
  85                         //initialize outputTypes
  86                         vector<string> tempOutNames;
  87                         outputTypes["fasta"] = tempOutNames;
  88                         outputTypes["accnos"] = tempOutNames;
  89
  90                         //if the user changes the input directory command factory will send this info to us in the output parameter
  91                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  92                         if (inputDir == "not found"){   inputDir = "";          }
  93                         else {
  94                                 string path;
  95                                 it = parameters.find("fasta");
  96                                 //user has given a template file
  97                                 if(it != parameters.end()){
  98                                         path = m->hasPath(it->second);
  99                                         //if the user has not given a path then, add inputdir. else leave path alone.
 100                                         if (path == "") {       parameters["fasta"] = inputDir + it->second;            }
 101                                 }
 102                         }
 103
 104                         //check for required parameters
 105                         fastafile = validParameter.validFile(parameters, "fasta", true);
 106                         if (fastafile == "not open") { abort = true; }
 107                         else if (fastafile == "not found") {  m->mothurOut("You must provide a fasta file."); m->mothurOutEndLine(); abort = true; }
 108
 109                         //if the user changes the output directory command factory will send this info to us in the output parameter
 110                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(fastafile);      }
 111
 112                         string temp = validParameter.validFile(parameters, "numbases", false);  if (temp == "not found") { temp = "0"; }
 113                         convert(temp, numbases);
 114
 115                         temp = validParameter.validFile(parameters, "countgaps", false);        if (temp == "not found") { temp = "f"; }
 116                         countGaps = m->isTrue(temp);
 117
 118                         temp = validParameter.validFile(parameters, "short", false);    if (temp == "not found") { temp = "f"; }
 119                         Short = m->isTrue(temp);
 120
 121                         keep = validParameter.validFile(parameters, "keep", false);             if (keep == "not found") { keep = "front"; }
 122
 123                         if (numbases == 0)  { m->mothurOut("You must provide the number of bases you want to keep for the chops.seqs command."); m->mothurOutEndLine(); abort = true;  }
 124                 }
 125
 126         }
 127         catch(exception& e) {
 128                 m->errorOut(e, "ChopSeqsCommand", "ChopSeqsCommand");
 129                 exit(1);
 130         }
 131 }
 132 //**********************************************************************************************************************
 133
 134 void ChopSeqsCommand::help(){
 135         try {
 136                 m->mothurOut("The chop.seqs command reads a fasta file and outputs a .chop.fasta containing the trimmed sequences. Note: If a sequence is completely 'chopped', an accnos file will be created with the names of the sequences removed. \n");
 137                 m->mothurOut("The chop.seqs command parameters are fasta, numbases, countgaps and keep. fasta and numbases are required required.\n");
 138                 m->mothurOut("The chop.seqs command should be in the following format: chop.seqs(fasta=yourFasta, numbases=yourNum, keep=yourKeep).\n");
 139                 m->mothurOut("The numbases parameter allows you to specify the number of bases you want to keep.\n");
 140                 m->mothurOut("The keep parameter allows you to specify whether you want to keep the front or the back of your sequence, default=front.\n");
 141                 m->mothurOut("The countgaps parameter allows you to specify whether you want to count gaps as bases, default=false.\n");
 142                 m->mothurOut("The short parameter allows you to specify you want to keep sequences that are too short to chop, default=false.\n");
 143                 m->mothurOut("For example, if you ran chop.seqs with numbases=200 and short=t, if a sequence had 100 bases mothur would keep the sequence rather than eliminate it.\n");
 144                 m->mothurOut("Example chop.seqs(fasta=amazon.fasta, numbases=200, keep=front).\n");
 145                 m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n");
 146         }
 147         catch(exception& e) {
 148                 m->errorOut(e, "ChopSeqsCommand", "help");
 149                 exit(1);
 150         }
 151 }
 152
 153 //**********************************************************************************************************************
 154
 155 int ChopSeqsCommand::execute(){
 156         try {
 157
 158                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 159
 160                 string outputFileName = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "chop.fasta";
 161                 string outputFileNameAccnos = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "chop.accnos";
 162
 163                 ofstream out;
 164                 m->openOutputFile(outputFileName, out);
 165
 166                 ofstream outAcc;
 167                 m->openOutputFile(outputFileNameAccnos, outAcc);
 168
 169                 ifstream in;
 170                 m->openInputFile(fastafile, in);
 171
 172                 bool wroteAccnos = false;
 173
 174                 while (!in.eof()) {
 175
 176                         Sequence seq(in);
 177
 178                         if (m->control_pressed) { outputTypes.clear(); in.close(); out.close(); outAcc.close(); remove(outputFileName.c_str()); remove(outputFileNameAccnos.c_str()); return 0;  }
 179
 180                         if (seq.getName() != "") {
 181                                 string newSeqString = getChopped(seq);
 182
 183                                 //output trimmed sequence
 184                                 if (newSeqString != "") {
 185                                         out << ">" << seq.getName() << endl << newSeqString << endl;
 186                                 }else{
 187                                         outAcc << seq.getName() << endl;
 188                                         wroteAccnos = true;
 189                                 }
 190                         }
 191                 }
 192                 in.close();
 193                 out.close();
 194                 outAcc.close();
 195
 196                 m->mothurOutEndLine();
 197                 m->mothurOut("Output File Name: "); m->mothurOutEndLine();
 198                 m->mothurOut(outputFileName); m->mothurOutEndLine();    outputNames.push_back(outputFileName); outputTypes["fasta"].push_back(outputFileName);
 199
 200                 if (wroteAccnos) { m->mothurOut(outputFileNameAccnos); m->mothurOutEndLine(); outputNames.push_back(outputFileNameAccnos); outputTypes["accnos"].push_back(outputFileNameAccnos); }
 201                 else {  remove(outputFileNameAccnos.c_str());  }
 202
 203                 m->mothurOutEndLine();
 204
 205                 //set fasta file as new current fastafile
 206                 string current = "";
 207                 itTypes = outputTypes.find("fasta");
 208                 if (itTypes != outputTypes.end()) {
 209                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
 210                 }
 211
 212                 if (wroteAccnos) { //set accnos file as new current accnosfile
 213                         itTypes = outputTypes.find("accnos");
 214                         if (itTypes != outputTypes.end()) {
 215                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setAccnosFile(current); }
 216                         }
 217                 }
 218
 219
 220                 return 0;
 221         }
 222
 223         catch(exception& e) {
 224                 m->errorOut(e, "ChopSeqsCommand", "execute");
 225                 exit(1);
 226         }
 227 }
 228 //**********************************************************************************************************************
 229 string ChopSeqsCommand::getChopped(Sequence seq) {
 230         try {
 231                 string temp = seq.getAligned();
 232                 string tempUnaligned = seq.getUnaligned();
 233
 234                 if (countGaps) {
 235                         //if needed trim sequence
 236                         if (keep == "front") {//you want to keep the beginning
 237                                 int tempLength = temp.length();
 238
 239                                 if (tempLength > numbases) { //you have enough bases to remove some
 240
 241                                         int stopSpot = 0;
 242                                         int numBasesCounted = 0;
 243
 244                                         for (int i = 0; i < temp.length(); i++) {
 245                                                 //eliminate N's
 246                                                 if (toupper(temp[i]) == 'N') { temp[i] == '.'; }
 247
 248                                                 numBasesCounted++;
 249
 250                                                 if (numBasesCounted >= numbases) { stopSpot = i; break; }
 251                                         }
 252
 253                                         if (stopSpot == 0) { temp = ""; }
 254                                         else {  temp = temp.substr(0, stopSpot);  }
 255
 256                                 }else {
 257                                         if (!Short) { temp = ""; } //sequence too short
 258                                 }
 259                         }else { //you are keeping the back
 260                                 int tempLength = temp.length();
 261                                 if (tempLength > numbases) { //you have enough bases to remove some
 262
 263                                         int stopSpot = 0;
 264                                         int numBasesCounted = 0;
 265
 266                                         for (int i = (temp.length()-1); i >= 0; i--) {
 267                                                 //eliminate N's
 268                                                 if (toupper(temp[i]) == 'N') { temp[i] == '.'; }
 269
 270                                                 numBasesCounted++;
 271
 272                                                 if (numBasesCounted >= numbases) { stopSpot = i; break; }
 273                                         }
 274
 275                                         if (stopSpot == 0) { temp = ""; }
 276                                         else {  temp = temp.substr(stopSpot+1);  }
 277                                 }else {
 278                                         if (!Short) { temp = ""; } //sequence too short
 279                                 }
 280                         }
 281
 282                 }else{
 283
 284                         //if needed trim sequence
 285                         if (keep == "front") {//you want to keep the beginning
 286                                 int tempLength = tempUnaligned.length();
 287
 288                                 if (tempLength > numbases) { //you have enough bases to remove some
 289
 290                                         int stopSpot = 0;
 291                                         int numBasesCounted = 0;
 292
 293                                         for (int i = 0; i < temp.length(); i++) {
 294                                                 //eliminate N's
 295                                                 if (toupper(temp[i]) == 'N') {
 296                                                         temp[i] == '.';
 297                                                         tempLength--;
 298                                                         if (tempLength < numbases) { stopSpot = 0; break; }
 299                                                 }
 300
 301                                                 if(isalpha(temp[i])) { numBasesCounted++; }
 302
 303                                                 if (numBasesCounted >= numbases) { stopSpot = i; break; }
 304                                         }
 305
 306                                         if (stopSpot == 0) { temp = ""; }
 307                                         else {  temp = temp.substr(0, stopSpot);  }
 308
 309                                 }else {
 310                                         if (!Short) { temp = ""; } //sequence too short
 311                                 }
 312                         }else { //you are keeping the back
 313                                 int tempLength = tempUnaligned.length();
 314                                 if (tempLength > numbases) { //you have enough bases to remove some
 315
 316                                         int stopSpot = 0;
 317                                         int numBasesCounted = 0;
 318
 319                                         for (int i = (temp.length()-1); i >= 0; i--) {
 320                                                 //eliminate N's
 321                                                 if (toupper(temp[i]) == 'N') {
 322                                                         temp[i] == '.';
 323                                                         tempLength--;
 324                                                         if (tempLength < numbases) { stopSpot = 0; break; }
 325                                                 }
 326
 327                                                 if(isalpha(temp[i])) { numBasesCounted++; }
 328
 329                                                 if (numBasesCounted >= numbases) { stopSpot = i; break; }
 330                                         }
 331
 332                                         if (stopSpot == 0) { temp = ""; }
 333                                         else {  temp = temp.substr(stopSpot+1);  }
 334                                 }else {
 335                                         if (!Short) { temp = ""; } //sequence too short
 336                                 }
 337                         }
 338                 }
 339
 340                 return temp;
 341         }
 342         catch(exception& e) {
 343                 m->errorOut(e, "ChopSeqsCommand", "getChopped");
 344                 exit(1);
 345         }
 346 }
 347 //**********************************************************************************************************************
 348
 349