screenseqscommand.cpp

   1 /*
   2  *  screenseqscommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Pat Schloss on 6/3/09.
   6  *  Copyright 2009 Patrick D. Schloss. All rights reserved.
   7  *
   8  */
   9
  10 #include "screenseqscommand.h"
  11 #include "sequence.hpp"
  12
  13 //***************************************************************************************************************
  14
  15 ScreenSeqsCommand::ScreenSeqsCommand(string option){
  16         try {
  17                 abort = false;
  18
  19                 //allow user to run help
  20                 if(option == "help") { help(); abort = true; }
  21
  22                 else {
  23                         //valid paramters for this command
  24                         string AlignArray[] =  {"fasta", "start", "end", "maxambig", "maxhomop", "minlength", "maxlength",
  25                                                                         "name", "group", "alignreport"};
  26                         vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
  27
  28                         OptionParser parser(option);
  29                         map<string,string> parameters = parser.getParameters();
  30
  31                         ValidParameters validParameter;
  32
  33                         //check to make sure all parameters are valid for command
  34                         for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
  35                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  36                         }
  37
  38                         //check for required parameters
  39                         fastafile = validParameter.validFile(parameters, "fasta", true);
  40                         if (fastafile == "not found") { mothurOut("fasta is a required parameter for the screen.seqs command."); mothurOutEndLine(); abort = true; }
  41                         else if (fastafile == "not open") { abort = true; }
  42
  43                         groupfile = validParameter.validFile(parameters, "group", true);
  44                         if (groupfile == "not open") { abort = true; }
  45                         else if (groupfile == "not found") { groupfile = ""; }
  46
  47                         namefile = validParameter.validFile(parameters, "name", true);
  48                         if (namefile == "not open") { abort = true; }
  49                         else if (namefile == "not found") { namefile = ""; }
  50
  51                         alignreport = validParameter.validFile(parameters, "alignreport", true);
  52                         if (alignreport == "not open") { abort = true; }
  53                         else if (alignreport == "not found") { alignreport = ""; }
  54
  55                         //check for optional parameter and set defaults
  56                         // ...at some point should added some additional type checking...
  57                         string temp;
  58                         temp = validParameter.validFile(parameters, "start", false);            if (temp == "not found") { temp = "-1"; }
  59                         convert(temp, startPos);
  60
  61                         temp = validParameter.validFile(parameters, "end", false);                      if (temp == "not found") { temp = "-1"; }
  62                         convert(temp, endPos);
  63
  64                         temp = validParameter.validFile(parameters, "maxambig", false);         if (temp == "not found") { temp = "-1"; }
  65                         convert(temp, maxAmbig);
  66
  67                         temp = validParameter.validFile(parameters, "maxhomop", false);         if (temp == "not found") { temp = "-1"; }
  68                         convert(temp, maxHomoP);
  69
  70                         temp = validParameter.validFile(parameters, "minlength", false);        if (temp == "not found") { temp = "-1"; }
  71                         convert(temp, minLength);
  72
  73                         temp = validParameter.validFile(parameters, "maxlength", false);        if (temp == "not found") { temp = "-1"; }
  74                         convert(temp, maxLength);
  75                 }
  76
  77         }
  78         catch(exception& e) {
  79                 errorOut(e, "ScreenSeqsCommand", "ScreenSeqsCommand");
  80                 exit(1);
  81         }
  82 }
  83 //**********************************************************************************************************************
  84
  85 void ScreenSeqsCommand::help(){
  86         try {
  87                 mothurOut("The screen.seqs command reads a fastafile and creates .....\n");
  88                 mothurOut("The screen.seqs command parameters are fasta, start, end, maxambig, maxhomop, minlength, maxlength, name, and group.\n");
  89                 mothurOut("The fasta parameter is required.\n");
  90                 mothurOut("The start parameter .... The default is -1.\n");
  91                 mothurOut("The end parameter .... The default is -1.\n");
  92                 mothurOut("The maxambig parameter .... The default is -1.\n");
  93                 mothurOut("The maxhomop parameter .... The default is -1.\n");
  94                 mothurOut("The minlength parameter .... The default is -1.\n");
  95                 mothurOut("The maxlength parameter .... The default is -1.\n");
  96                 mothurOut("The name parameter allows you to provide a namesfile, and the group parameter allows you to provide a groupfile.\n");
  97                 mothurOut("The screen.seqs command should be in the following format: \n");
  98                 mothurOut("screen.seqs(fasta=yourFastaFile, name=youNameFile, group=yourGroupFIle, start=yourStart, end=yourEnd, maxambig=yourMaxambig,  \n");
  99                 mothurOut("maxhomop=yourMaxhomop, minlength=youMinlength, maxlength=yourMaxlength)  \n");
 100                 mothurOut("Example screen.seqs(fasta=abrecovery.fasta, name=abrecovery.names, group=abrecovery.groups, start=..., end=..., maxambig=..., maxhomop=..., minlength=..., maxlength=...).\n");
 101                 mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n");
 102
 103         }
 104         catch(exception& e) {
 105                 errorOut(e, "ScreenSeqsCommand", "help");
 106                 exit(1);
 107         }
 108 }
 109
 110 //***************************************************************************************************************
 111
 112 ScreenSeqsCommand::~ScreenSeqsCommand(){        /*      do nothing      */      }
 113
 114 //***************************************************************************************************************
 115
 116 int ScreenSeqsCommand::execute(){
 117         try{
 118
 119                 if (abort == true) { return 0; }
 120
 121                 ifstream inFASTA;
 122                 openInputFile(fastafile, inFASTA);
 123
 124                 set<string> badSeqNames;
 125
 126                 string goodSeqFile = getRootName(fastafile) + "good" + getExtension(fastafile);
 127                 string badSeqFile = getRootName(fastafile) + "bad" + getExtension(fastafile);
 128
 129                 ofstream goodSeqOut;    openOutputFile(goodSeqFile, goodSeqOut);
 130                 ofstream badSeqOut;             openOutputFile(badSeqFile, badSeqOut);
 131
 132                 while(!inFASTA.eof()){
 133                         Sequence currSeq(inFASTA);
 134                         if (currSeq.getName() != "") {
 135                                 bool goodSeq = 1;               //      innocent until proven guilty
 136                                 if(goodSeq == 1 && startPos != -1 && startPos < currSeq.getStartPos())                  {       goodSeq = 0;    }
 137                                 if(goodSeq == 1 && endPos != -1 && endPos > currSeq.getEndPos())                                {       goodSeq = 0;    }
 138                                 if(goodSeq == 1 && maxAmbig != -1 && maxAmbig < currSeq.getAmbigBases())                {       goodSeq = 0;    }
 139                                 if(goodSeq == 1 && maxHomoP != -1 && maxHomoP < currSeq.getLongHomoPolymer())   {       goodSeq = 0;    }
 140                                 if(goodSeq == 1 && minLength != -1 && minLength > currSeq.getNumBases())                {       goodSeq = 0;    }
 141                                 if(goodSeq == 1 && maxLength != -1 && maxLength < currSeq.getNumBases())                {       goodSeq = 0;    }
 142
 143                                 if(goodSeq == 1){
 144                                         currSeq.printSequence(goodSeqOut);
 145                                 }
 146                                 else{
 147                                         currSeq.printSequence(badSeqOut);
 148                                         badSeqNames.insert(currSeq.getName());
 149                                 }
 150                         }
 151                         gobble(inFASTA);
 152                 }
 153                 if(namefile != "" && groupfile != "")   {       screenNameGroupFile(badSeqNames);       }       // this screens both names and groups
 154                 else if(namefile != "")                                 {       screenNameGroupFile(badSeqNames);       }
 155                 else if(groupfile != "")                                {       screenGroupFile(badSeqNames);           }       // this screens just the groups
 156                 if(alignreport != "")                                   {       screenAlignReport(badSeqNames);         }
 157
 158                 goodSeqOut.close();
 159                 badSeqOut.close();
 160                 inFASTA.close();
 161                 return 0;
 162         }
 163         catch(exception& e) {
 164                 errorOut(e, "ScreenSeqsCommand", "execute");
 165                 exit(1);
 166         }
 167 }
 168
 169 //***************************************************************************************************************
 170
 171 void ScreenSeqsCommand::screenNameGroupFile(set<string> badSeqNames){
 172
 173         ifstream inputNames;
 174         openInputFile(namefile, inputNames);
 175         set<string> badSeqGroups;
 176         string seqName, seqList, group;
 177         set<string>::iterator it;
 178
 179         string goodNameFile = getRootName(namefile) + "good" + getExtension(namefile);
 180         string badNameFile = getRootName(namefile) + "bad" + getExtension(namefile);
 181
 182         ofstream goodNameOut;   openOutputFile(goodNameFile, goodNameOut);
 183         ofstream badNameOut;    openOutputFile(badNameFile, badNameOut);
 184
 185         while(!inputNames.eof()){
 186                 inputNames >> seqName >> seqList;
 187                 it = badSeqNames.find(seqName);
 188
 189                 if(it != badSeqNames.end()){
 190                         badSeqNames.erase(it);
 191                         badNameOut << seqName << '\t' << seqList << endl;
 192                         if(namefile != ""){
 193                                 int start = 0;
 194                                 for(int i=0;i<seqList.length();i++){
 195                                         if(seqList[i] == ','){
 196                                                 badSeqGroups.insert(seqList.substr(start,i-start));
 197                                                 start = i+1;
 198                                         }
 199                                 }
 200                                 badSeqGroups.insert(seqList.substr(start,seqList.length()-start));
 201                         }
 202                 }
 203                 else{
 204                         goodNameOut << seqName << '\t' << seqList << endl;
 205                 }
 206                 gobble(inputNames);
 207         }
 208         inputNames.close();
 209         goodNameOut.close();
 210         badNameOut.close();
 211
 212         //we were unable to remove some of the bad sequences
 213         if (badSeqNames.size() != 0) {
 214                 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
 215                         mothurOut("Your namefile does not include the sequence " + *it + " please correct.");
 216                         mothurOutEndLine();
 217                 }
 218         }
 219
 220         if(groupfile != ""){
 221
 222                 ifstream inputGroups;
 223                 openInputFile(groupfile, inputGroups);
 224
 225                 string goodGroupFile = getRootName(groupfile) + "good" + getExtension(groupfile);
 226                 string badGroupFile = getRootName(groupfile) + "bad" + getExtension(groupfile);
 227
 228                 ofstream goodGroupOut;  openOutputFile(goodGroupFile, goodGroupOut);
 229                 ofstream badGroupOut;   openOutputFile(badGroupFile, badGroupOut);
 230
 231                 while(!inputGroups.eof()){
 232                         inputGroups >> seqName >> group;
 233
 234                         it = badSeqGroups.find(seqName);
 235
 236                         if(it != badSeqGroups.end()){
 237                                 badSeqGroups.erase(it);
 238                                 badGroupOut << seqName << '\t' << group << endl;
 239                         }
 240                         else{
 241                                 goodGroupOut << seqName << '\t' << group << endl;
 242                         }
 243                         gobble(inputGroups);
 244                 }
 245                 inputGroups.close();
 246                 goodGroupOut.close();
 247                 badGroupOut.close();
 248
 249                 //we were unable to remove some of the bad sequences
 250                 if (badSeqGroups.size() != 0) {
 251                         for (it = badSeqGroups.begin(); it != badSeqGroups.end(); it++) {
 252                                 mothurOut("Your namefile does not include the sequence " + *it + " please correct.");
 253                                 mothurOutEndLine();
 254                         }
 255                 }
 256         }
 257 }
 258
 259 //***************************************************************************************************************
 260
 261 void ScreenSeqsCommand::screenGroupFile(set<string> badSeqNames){
 262
 263         ifstream inputGroups;
 264         openInputFile(groupfile, inputGroups);
 265         string seqName, group;
 266         set<string>::iterator it;
 267
 268         string goodGroupFile = getRootName(groupfile) + "good" + getExtension(groupfile);
 269         string badGroupFile = getRootName(groupfile) + "bad" + getExtension(groupfile);
 270
 271         ofstream goodGroupOut;  openOutputFile(goodGroupFile, goodGroupOut);
 272         ofstream badGroupOut;   openOutputFile(badGroupFile, badGroupOut);
 273
 274         while(!inputGroups.eof()){
 275                 inputGroups >> seqName >> group;
 276                 it = badSeqNames.find(seqName);
 277
 278                 if(it != badSeqNames.end()){
 279                         badSeqNames.erase(it);
 280                         badGroupOut << seqName << '\t' << group << endl;
 281                 }
 282                 else{
 283                         goodGroupOut << seqName << '\t' << group << endl;
 284                 }
 285                 gobble(inputGroups);
 286         }
 287
 288         //we were unable to remove some of the bad sequences
 289         if (badSeqNames.size() != 0) {
 290                 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
 291                         mothurOut("Your groupfile does not include the sequence " + *it + " please correct.");
 292                         mothurOutEndLine();
 293                 }
 294         }
 295
 296         inputGroups.close();
 297         goodGroupOut.close();
 298         badGroupOut.close();
 299
 300 }
 301
 302 //***************************************************************************************************************
 303
 304 void ScreenSeqsCommand::screenAlignReport(set<string> badSeqNames){
 305
 306         ifstream inputAlignReport;
 307         openInputFile(alignreport, inputAlignReport);
 308         string seqName, group;
 309         set<string>::iterator it;
 310
 311         string goodAlignReportFile = getRootName(alignreport) + "good" + getExtension(alignreport);
 312         string badAlignReportFile = getRootName(alignreport) + "bad" + getExtension(alignreport);
 313
 314         ofstream goodAlignReportOut;    openOutputFile(goodAlignReportFile, goodAlignReportOut);
 315         ofstream badAlignReportOut;             openOutputFile(badAlignReportFile, badAlignReportOut);
 316
 317         while (!inputAlignReport.eof()) {               //      need to copy header
 318                 char c = inputAlignReport.get();
 319                 goodAlignReportOut << c;
 320                 badAlignReportOut << c;
 321                 if (c == 10 || c == 13){        break;  }
 322         }
 323
 324         while(!inputAlignReport.eof()){
 325                 inputAlignReport >> seqName;
 326                 it = badSeqNames.find(seqName);
 327                 string line;
 328                 while (!inputAlignReport.eof()) {               //      need to copy header
 329                         char c = inputAlignReport.get();
 330                         line += c;
 331                         if (c == 10 || c == 13){        break;  }
 332                 }
 333
 334                 if(it != badSeqNames.end()){
 335                         badSeqNames.erase(it);
 336                         badAlignReportOut << seqName << '\t' << line;;
 337                 }
 338                 else{
 339                         goodAlignReportOut << seqName << '\t' << line;
 340                 }
 341                 gobble(inputAlignReport);
 342         }
 343
 344         //we were unable to remove some of the bad sequences
 345         if (badSeqNames.size() != 0) {
 346                 for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
 347                         mothurOut("Your file does not include the sequence " + *it + " please correct.");
 348                         mothurOutEndLine();
 349                 }
 350         }
 351
 352         inputAlignReport.close();
 353         goodAlignReportOut.close();
 354         badAlignReportOut.close();
 355
 356 }
 357
 358 //***************************************************************************************************************
 359
 360