]> git.donarmstrong.com Git - mothur.git/blobdiff - screenseqscommand.cpp
pat's mods to morisitahorn and pre.cluster
[mothur.git] / screenseqscommand.cpp
index 5789fd0665f849822ec8d8c365de3a26e65b77ef..32426f3092d17468aa2f2599528f915423954376 100644 (file)
@@ -21,7 +21,8 @@ ScreenSeqsCommand::ScreenSeqsCommand(string option){
                
                else {
                        //valid paramters for this command
-                       string AlignArray[] =  {"fasta", "start", "end", "maxambig", "maxhomop", "minlength", "maxlength", "name", "group"};
+                       string AlignArray[] =  {"fasta", "start", "end", "maxambig", "maxhomop", "minlength", "maxlength",
+                                                                       "name", "group", "alignreport"};
                        vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
                        
                        OptionParser parser(option);
@@ -36,9 +37,9 @@ ScreenSeqsCommand::ScreenSeqsCommand(string option){
                        
                        //check for required parameters
                        fastafile = validParameter.validFile(parameters, "fasta", true);
-                       if (fastafile == "not found") { cout << "fasta is a required parameter for the screen.seqs command." << endl; abort = true; }
+                       if (fastafile == "not found") { mothurOut("fasta is a required parameter for the screen.seqs command."); mothurOutEndLine(); abort = true; }
                        else if (fastafile == "not open") { abort = true; }     
-               
+       
                        groupfile = validParameter.validFile(parameters, "group", true);
                        if (groupfile == "not open") { abort = true; }  
                        else if (groupfile == "not found") { groupfile = ""; }
@@ -47,14 +48,17 @@ ScreenSeqsCommand::ScreenSeqsCommand(string option){
                        if (namefile == "not open") { abort = true; }
                        else if (namefile == "not found") { namefile = ""; }    
 
-               
+                       alignreport = validParameter.validFile(parameters, "alignreport", true);
+                       if (alignreport == "not open") { abort = true; }
+                       else if (alignreport == "not found") { alignreport = ""; }      
+                       
                        //check for optional parameter and set defaults
                        // ...at some point should added some additional type checking...
                        string temp;
-                       temp = validParameter.validFile(parameters, "start", false);                    if (temp == "not found") { temp = "-1"; }
+                       temp = validParameter.validFile(parameters, "start", false);            if (temp == "not found") { temp = "-1"; }
                        convert(temp, startPos); 
                
-                       temp = validParameter.validFile(parameters, "end", false);                              if (temp == "not found") { temp = "-1"; }
+                       temp = validParameter.validFile(parameters, "end", false);                      if (temp == "not found") { temp = "-1"; }
                        convert(temp, endPos);  
 
                        temp = validParameter.validFile(parameters, "maxambig", false);         if (temp == "not found") { temp = "-1"; }
@@ -63,52 +67,44 @@ ScreenSeqsCommand::ScreenSeqsCommand(string option){
                        temp = validParameter.validFile(parameters, "maxhomop", false);         if (temp == "not found") { temp = "-1"; }
                        convert(temp, maxHomoP);  
 
-                       temp = validParameter.validFile(parameters, "minlength", false);                if (temp == "not found") { temp = "-1"; }
+                       temp = validParameter.validFile(parameters, "minlength", false);        if (temp == "not found") { temp = "-1"; }
                        convert(temp, minLength); 
                        
-                       temp = validParameter.validFile(parameters, "maxlength", false);                if (temp == "not found") { temp = "-1"; }
+                       temp = validParameter.validFile(parameters, "maxlength", false);        if (temp == "not found") { temp = "-1"; }
                        convert(temp, maxLength); 
                }
 
        }
        catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the ScreenSeqsCommand class Function ScreenSeqsCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               errorOut(e, "ScreenSeqsCommand", "ScreenSeqsCommand");
                exit(1);
        }
-       catch(...) {
-               cout << "An unknown error has occurred in the ScreenSeqsCommand class function ScreenSeqsCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }       
 }
 //**********************************************************************************************************************
 
 void ScreenSeqsCommand::help(){
        try {
-               cout << "The screen.seqs command reads a fastafile and creates ....." << "\n";
-               cout << "The screen.seqs command parameters are fasta, start, end, maxambig, maxhomop, minlength, maxlength, name, and group." << "\n";
-               cout << "The fasta parameter is required." << "\n";
-               cout << "The start parameter .... The default is -1." << "\n";
-               cout << "The end parameter .... The default is -1." << "\n";
-               cout << "The maxambig parameter .... The default is -1." << "\n";
-               cout << "The maxhomop parameter .... The default is -1." << "\n";
-               cout << "The minlength parameter .... The default is -1." << "\n";
-               cout << "The maxlength parameter .... The default is -1." << "\n";
-               cout << "The name parameter allows you to provide a namesfile, and the group parameter allows you to provide a groupfile." << "\n";
-               cout << "The screen.seqs command should be in the following format: " << "\n";
-               cout << "screen.seqs(fasta=yourFastaFile, name=youNameFile, group=yourGroupFIle, start=yourStart, end=yourEnd, maxambig=yourMaxambig,  " << "\n";
-               cout << "maxhomop=yourMaxhomop, minlength=youMinlength, maxlength=yourMaxlength)  " << "\n";    
-               cout << "Example screen.seqs(fasta=abrecovery.fasta, name=abrecovery.names, group=abrecovery.groups, start=..., end=..., maxambig=..., maxhomop=..., minlength=..., maxlength=...)." << "\n";
-               cout << "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta)." << "\n" << "\n";
+               mothurOut("The screen.seqs command reads a fastafile and creates .....\n");
+               mothurOut("The screen.seqs command parameters are fasta, start, end, maxambig, maxhomop, minlength, maxlength, name, and group.\n");
+               mothurOut("The fasta parameter is required.\n");
+               mothurOut("The start parameter .... The default is -1.\n");
+               mothurOut("The end parameter .... The default is -1.\n");
+               mothurOut("The maxambig parameter .... The default is -1.\n");
+               mothurOut("The maxhomop parameter .... The default is -1.\n");
+               mothurOut("The minlength parameter .... The default is -1.\n");
+               mothurOut("The maxlength parameter .... The default is -1.\n");
+               mothurOut("The name parameter allows you to provide a namesfile, and the group parameter allows you to provide a groupfile.\n");
+               mothurOut("The screen.seqs command should be in the following format: \n");
+               mothurOut("screen.seqs(fasta=yourFastaFile, name=youNameFile, group=yourGroupFIle, start=yourStart, end=yourEnd, maxambig=yourMaxambig,  \n");
+               mothurOut("maxhomop=yourMaxhomop, minlength=youMinlength, maxlength=yourMaxlength)  \n");       
+               mothurOut("Example screen.seqs(fasta=abrecovery.fasta, name=abrecovery.names, group=abrecovery.groups, start=..., end=..., maxambig=..., maxhomop=..., minlength=..., maxlength=...).\n");
+               mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n");
 
        }
        catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the ScreenSeqsCommand class Function help. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               errorOut(e, "ScreenSeqsCommand", "help");
                exit(1);
        }
-       catch(...) {
-               cout << "An unknown error has occurred in the ScreenSeqsCommand class function help. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }       
 }
 
 //***************************************************************************************************************
@@ -135,41 +131,39 @@ int ScreenSeqsCommand::execute(){
                
                while(!inFASTA.eof()){
                        Sequence currSeq(inFASTA);
-                       bool goodSeq = 1;               //      innocent until proven guilty
-                       if(goodSeq == 1 && startPos != -1 && startPos < currSeq.getStartPos())                  {       goodSeq = 0;    }
-                       if(goodSeq == 1 && endPos != -1 && endPos > currSeq.getEndPos())                                {       goodSeq = 0;    }
-                       if(goodSeq == 1 && maxAmbig != -1 && maxAmbig < currSeq.getAmbigBases())                {       goodSeq = 0;    }
-                       if(goodSeq == 1 && maxHomoP != -1 && maxHomoP < currSeq.getLongHomoPolymer())   {       goodSeq = 0;    }
-                       if(goodSeq == 1 && minLength != -1 && minLength > currSeq.getNumBases())                {       goodSeq = 0;    }
-                       if(goodSeq == 1 && maxLength != -1 && maxLength < currSeq.getNumBases())                {       goodSeq = 0;    }
-                       
-                       if(goodSeq == 1){
-                               currSeq.printSequence(goodSeqOut);      
-                       }
-                       else{
-                               currSeq.printSequence(badSeqOut);       
-                               badSeqNames.insert(currSeq.getName());
+                       if (currSeq.getName() != "") {
+                               bool goodSeq = 1;               //      innocent until proven guilty
+                               if(goodSeq == 1 && startPos != -1 && startPos < currSeq.getStartPos())                  {       goodSeq = 0;    }
+                               if(goodSeq == 1 && endPos != -1 && endPos > currSeq.getEndPos())                                {       goodSeq = 0;    }
+                               if(goodSeq == 1 && maxAmbig != -1 && maxAmbig < currSeq.getAmbigBases())                {       goodSeq = 0;    }
+                               if(goodSeq == 1 && maxHomoP != -1 && maxHomoP < currSeq.getLongHomoPolymer())   {       goodSeq = 0;    }
+                               if(goodSeq == 1 && minLength != -1 && minLength > currSeq.getNumBases())                {       goodSeq = 0;    }
+                               if(goodSeq == 1 && maxLength != -1 && maxLength < currSeq.getNumBases())                {       goodSeq = 0;    }
+                               
+                               if(goodSeq == 1){
+                                       currSeq.printSequence(goodSeqOut);      
+                               }
+                               else{
+                                       currSeq.printSequence(badSeqOut);       
+                                       badSeqNames.insert(currSeq.getName());
+                               }
                        }
                        gobble(inFASTA);
                }       
-               if(namefile != ""){
-                       screenNameGroupFile(badSeqNames);
-               }
-               else if(groupfile != ""){
-                       screenGroupFile(badSeqNames);
-               }
+               if(namefile != "" && groupfile != "")   {       screenNameGroupFile(badSeqNames);       }       // this screens both names and groups
+               else if(namefile != "")                                 {       screenNameGroupFile(badSeqNames);       }
+               else if(groupfile != "")                                {       screenGroupFile(badSeqNames);           }       // this screens just the groups
+               if(alignreport != "")                                   {       screenAlignReport(badSeqNames);         }
                
+               goodSeqOut.close();
+               badSeqOut.close();
+               inFASTA.close();
                return 0;
        }
        catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the ScreenSeqsCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               errorOut(e, "ScreenSeqsCommand", "execute");
                exit(1);
        }
-       catch(...) {
-               cout << "An unknown error has occurred in the ScreenSeqsCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }
-       
 }
 
 //***************************************************************************************************************
@@ -215,6 +209,14 @@ void ScreenSeqsCommand::screenNameGroupFile(set<string> badSeqNames){
        goodNameOut.close();
        badNameOut.close();
        
+       //we were unable to remove some of the bad sequences
+       if (badSeqNames.size() != 0) {
+               for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {  
+                       mothurOut("Your namefile does not include the sequence " + *it + " please correct."); 
+                       mothurOutEndLine();
+               }
+       }
+
        if(groupfile != ""){
                
                ifstream inputGroups;
@@ -243,6 +245,14 @@ void ScreenSeqsCommand::screenNameGroupFile(set<string> badSeqNames){
                inputGroups.close();
                goodGroupOut.close();
                badGroupOut.close();
+               
+               //we were unable to remove some of the bad sequences
+               if (badSeqGroups.size() != 0) {
+                       for (it = badSeqGroups.begin(); it != badSeqGroups.end(); it++) {  
+                               mothurOut("Your namefile does not include the sequence " + *it + " please correct."); 
+                               mothurOutEndLine();
+                       }
+               }
        }
 }
 
@@ -274,6 +284,15 @@ void ScreenSeqsCommand::screenGroupFile(set<string> badSeqNames){
                }
                gobble(inputGroups);
        }
+       
+       //we were unable to remove some of the bad sequences
+       if (badSeqNames.size() != 0) {
+               for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {  
+                       mothurOut("Your groupfile does not include the sequence " + *it + " please correct."); 
+                       mothurOutEndLine();
+               }
+       }
+       
        inputGroups.close();
        goodGroupOut.close();
        badGroupOut.close();
@@ -282,4 +301,60 @@ void ScreenSeqsCommand::screenGroupFile(set<string> badSeqNames){
 
 //***************************************************************************************************************
 
+void ScreenSeqsCommand::screenAlignReport(set<string> badSeqNames){
+       
+       ifstream inputAlignReport;
+       openInputFile(alignreport, inputAlignReport);
+       string seqName, group;
+       set<string>::iterator it;
+       
+       string goodAlignReportFile = getRootName(alignreport) + "good" + getExtension(alignreport);
+       string badAlignReportFile = getRootName(alignreport) + "bad" + getExtension(alignreport);
+       
+       ofstream goodAlignReportOut;    openOutputFile(goodAlignReportFile, goodAlignReportOut);
+       ofstream badAlignReportOut;             openOutputFile(badAlignReportFile, badAlignReportOut);          
+
+       while (!inputAlignReport.eof()) {               //      need to copy header
+               char c = inputAlignReport.get();
+               goodAlignReportOut << c;
+               badAlignReportOut << c;
+               if (c == 10 || c == 13){        break;  }       
+       }
+
+       while(!inputAlignReport.eof()){
+               inputAlignReport >> seqName;
+               it = badSeqNames.find(seqName);
+               string line;            
+               while (!inputAlignReport.eof()) {               //      need to copy header
+                       char c = inputAlignReport.get();
+                       line += c;
+                       if (c == 10 || c == 13){        break;  }       
+               }
+               
+               if(it != badSeqNames.end()){
+                       badSeqNames.erase(it);
+                       badAlignReportOut << seqName << '\t' << line;;
+               }
+               else{
+                       goodAlignReportOut << seqName << '\t' << line;
+               }
+               gobble(inputAlignReport);
+       }
+       
+       //we were unable to remove some of the bad sequences
+       if (badSeqNames.size() != 0) {
+               for (it = badSeqNames.begin(); it != badSeqNames.end(); it++) {  
+                       mothurOut("Your file does not include the sequence " + *it + " please correct."); 
+                       mothurOutEndLine();
+               }
+       }
+
+       inputAlignReport.close();
+       goodAlignReportOut.close();
+       badAlignReportOut.close();
+       
+}
+
+//***************************************************************************************************************
+