]> git.donarmstrong.com Git - mothur.git/blobdiff - prcseqscommand.cpp
sffinfo bug with flow grams right index when clipQualRight=0
[mothur.git] / prcseqscommand.cpp
index d9c37769b07d24f40efbfef9e076da33f5c86c29..cf93cb796104aa90428e5466546c9839789216f7 100644 (file)
 //**********************************************************************************************************************
 vector<string> PcrSeqsCommand::setParameters(){        
        try {
-               CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pfasta);
-               CommandParameter poligos("oligos", "InputTypes", "", "", "ecolioligos", "none", "none",false,false); parameters.push_back(poligos);
-               CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
-        CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pgroup);
-        CommandParameter ptax("taxonomy", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(ptax);
-        CommandParameter pecoli("ecoli", "InputTypes", "", "", "ecolioligos", "none", "none",false,false); parameters.push_back(pecoli);
-               CommandParameter pstart("start", "Number", "", "-1", "", "", "",false,false); parameters.push_back(pstart);
-               CommandParameter pend("end", "Number", "", "-1", "", "", "",false,false); parameters.push_back(pend);
-               CommandParameter pnomatch("nomatch", "Multiple", "reject-keep", "reject", "", "", "",false,false); parameters.push_back(pnomatch);
-               CommandParameter ppdiffs("pdiffs", "Number", "", "0", "", "", "",false,false); parameters.push_back(ppdiffs);
-               CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
-               CommandParameter pkeepprimer("keepprimer", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pkeepprimer);
-        CommandParameter pkeepdots("keepdots", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pkeepdots);
-        CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
-               CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
+               CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none","fasta",false,true,true); parameters.push_back(pfasta);
+               CommandParameter poligos("oligos", "InputTypes", "", "", "ecolioligos", "none", "none","",false,false,true); parameters.push_back(poligos);
+        CommandParameter pname("name", "InputTypes", "", "", "NameCount", "none", "none","name",false,false,true); parameters.push_back(pname);
+        CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "none", "none","count",false,false,true); parameters.push_back(pcount);
+               CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "none", "none","group",false,false,true); parameters.push_back(pgroup);
+        CommandParameter ptax("taxonomy", "InputTypes", "", "", "none", "none", "none","taxonomy",false,false,true); parameters.push_back(ptax);
+        CommandParameter pecoli("ecoli", "InputTypes", "", "", "ecolioligos", "none", "none","",false,false); parameters.push_back(pecoli);
+               CommandParameter pstart("start", "Number", "", "-1", "", "", "","",false,false); parameters.push_back(pstart);
+               CommandParameter pend("end", "Number", "", "-1", "", "", "","",false,false); parameters.push_back(pend);
+               CommandParameter pnomatch("nomatch", "Multiple", "reject-keep", "reject", "", "", "","",false,false); parameters.push_back(pnomatch);
+        CommandParameter ppdiffs("pdiffs", "Number", "", "0", "", "", "","",false,false,true); parameters.push_back(ppdiffs);
+
+               CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
+               CommandParameter pkeepprimer("keepprimer", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(pkeepprimer);
+        CommandParameter pkeepdots("keepdots", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(pkeepdots);
+        CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
+               CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
         
                vector<string> myArray;
                for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
@@ -40,8 +42,16 @@ vector<string> PcrSeqsCommand::setParameters(){
 string PcrSeqsCommand::getHelpString(){        
        try {
                string helpString = "";
-               helpString += "The pcr.seqs command reads a fasta file ...\n";
-               
+               helpString += "The pcr.seqs command reads a fasta file.\n";
+        helpString += "The pcr.seqs command parameters are fasta, oligos, name, group, count, taxonomy, ecoli, start, end, nomatch, pdiffs, processors, keepprimer and keepdots.\n";
+               helpString += "The ecoli parameter is used to provide a fasta file containing a single reference sequence (e.g. for e. coli) this must be aligned. Mothur will trim to the start and end positions of the reference sequence.\n";
+        helpString += "The start parameter allows you to provide a starting position to trim to.\n";
+        helpString += "The end parameter allows you to provide a ending position to trim from.\n";
+        helpString += "The nomatch parameter allows you to decide what to do with sequences where the primer is not found. Default=reject, meaning remove from fasta file.  if nomatch=true, then do nothing to sequence.\n";
+        helpString += "The processors parameter allows you to use multiple processors.\n";
+        helpString += "The keepprimer parameter allows you to keep the primer, default=false.\n";
+        helpString += "The keepdots parameter allows you to keep the leading and trailing .'s, default=true.\n";
+        helpString += "The pdiffs parameter is used to specify the number of differences allowed in the primer. The default is 0.\n";
                helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n";
                helpString += "For more details please check out the wiki http://www.mothur.org/wiki/Pcr.seqs .\n";
                return helpString;
@@ -51,8 +61,26 @@ string PcrSeqsCommand::getHelpString(){
                exit(1);
        }
 }
-
-
+//**********************************************************************************************************************
+string PcrSeqsCommand::getOutputPattern(string type) {
+    try {
+        string pattern = "";
+        
+        if (type == "fasta")            {   pattern = "[filename],pcr,[extension]-[filename],[tag],pcr,[extension]";    }
+        else if (type == "taxonomy")    {   pattern = "[filename],pcr,[extension]";    }
+        else if (type == "name")        {   pattern = "[filename],pcr,[extension]";    }
+        else if (type == "group")       {   pattern = "[filename],pcr,[extension]";    }
+        else if (type == "count")       {   pattern = "[filename],pcr,[extension]";    }
+        else if (type == "accnos")      {   pattern = "[filename],bad.accnos";    }
+        else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
+        
+        return pattern;
+    }
+    catch(exception& e) {
+        m->errorOut(e, "PcrSeqsCommand", "getOutputPattern");
+        exit(1);
+    }
+}
 //**********************************************************************************************************************
 
 PcrSeqsCommand::PcrSeqsCommand(){      
@@ -64,6 +92,7 @@ PcrSeqsCommand::PcrSeqsCommand(){
                outputTypes["taxonomy"] = tempOutNames;
                outputTypes["group"] = tempOutNames;
                outputTypes["name"] = tempOutNames;
+        outputTypes["count"] = tempOutNames;
         outputTypes["accnos"] = tempOutNames;
        }
        catch(exception& e) {
@@ -103,6 +132,7 @@ PcrSeqsCommand::PcrSeqsCommand(string option)  {
                        outputTypes["group"] = tempOutNames;
                        outputTypes["name"] = tempOutNames;
             outputTypes["accnos"] = tempOutNames;
+            outputTypes["count"] = tempOutNames;
                        
                        //if the user changes the input directory command factory will send this info to us in the output parameter 
                        string inputDir = validParameter.validFile(parameters, "inputdir", false);              
@@ -156,11 +186,17 @@ PcrSeqsCommand::PcrSeqsCommand(string option)  {
                                        //if the user has not given a path then, add inputdir. else leave path alone.
                                        if (path == "") {       parameters["group"] = inputDir + it->second;            }
                                }
+                
+                it = parameters.find("count");
+                               //user has given a template file
+                               if(it != parameters.end()){ 
+                                       path = m->hasPath(it->second);
+                                       //if the user has not given a path then, add inputdir. else leave path alone.
+                                       if (path == "") {       parameters["count"] = inputDir + it->second;            }
+                               }
                                
                        }
             
-            //if the user changes the output directory command factory will send this info to us in the output parameter 
-                       outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = ""; }
                        
                        //check for required parameters
                        fastafile = validParameter.validFile(parameters, "fasta", true);
@@ -171,7 +207,9 @@ PcrSeqsCommand::PcrSeqsCommand(string option)  {
                        }else if (fastafile == "not open") { fastafile = ""; abort = true; }    
                        else { m->setFastaFile(fastafile); }
                        
-            
+            //if the user changes the output directory command factory will send this info to us in the output parameter 
+                       outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(fastafile);      }
+
                        //check for optional parameter and set defaults
                        // ...at some point should added some additional type checking...
                        string temp;
@@ -200,14 +238,24 @@ PcrSeqsCommand::PcrSeqsCommand(string option)  {
                        else if(groupfile == "not open"){       groupfile = ""; abort = true;   } 
             else { m->setGroupFile(groupfile); }
             
+            countfile = validParameter.validFile(parameters, "count", true);
+                       if (countfile == "not open") { countfile = ""; abort = true; }
+                       else if (countfile == "not found") { countfile = "";  } 
+                       else { m->setCountTableFile(countfile); }
+            
+            if ((namefile != "") && (countfile != "")) {
+                m->mothurOut("[ERROR]: you may only use one of the following: name or count."); m->mothurOutEndLine(); abort = true;
+            }
+                       
+            if ((groupfile != "") && (countfile != "")) {
+                m->mothurOut("[ERROR]: you may only use one of the following: group or count."); m->mothurOutEndLine(); abort=true;
+            }
+            
             taxfile = validParameter.validFile(parameters, "taxonomy", true);
                        if (taxfile == "not found"){    taxfile = "";           }
                        else if(taxfile == "not open"){ taxfile = ""; abort = true;     } 
             else { m->setTaxonomyFile(taxfile); }
-                       
-                       temp = validParameter.validFile(parameters, "pdiffs", false);           if (temp == "not found") { temp = "0"; }
-                       m->mothurConvert(temp, pdiffs);
-                       
+                                               
                        temp = validParameter.validFile(parameters, "start", false);    if (temp == "not found") { temp = "-1"; }
                        m->mothurConvert(temp, start);
             
@@ -217,6 +265,9 @@ PcrSeqsCommand::PcrSeqsCommand(string option)  {
                        temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = m->getProcessors();      }
                        m->setProcessors(temp);
                        m->mothurConvert(temp, processors); 
+            
+            temp = validParameter.validFile(parameters, "pdiffs", false);              if (temp == "not found") { temp = "0"; }
+                       m->mothurConvert(temp, pdiffs);
                        
             nomatch = validParameter.validFile(parameters, "nomatch", false);  if (nomatch == "not found") { nomatch = "reject"; }
                        
@@ -239,10 +290,12 @@ PcrSeqsCommand::PcrSeqsCommand(string option)  {
             }
                        
                        //check to make sure you didn't forget the name file by mistake                 
-                       if (namefile == "") {
-                               vector<string> files; files.push_back(fastafile);
-                               parser.getNameFile(files);
-                       }
+                       if (countfile == "") { 
+                if (namefile == "") {
+                    vector<string> files; files.push_back(fastafile);
+                    parser.getNameFile(files);
+                }
+            }
                }
         
        }
@@ -262,14 +315,17 @@ int PcrSeqsCommand::execute(){
         
         string thisOutputDir = outputDir;
                if (outputDir == "") {  thisOutputDir += m->hasPath(fastafile);  }
-               string trimSeqFile = thisOutputDir + m->getRootName(m->getSimpleName(fastafile)) + "pcr.fasta";
+        map<string, string> variables; 
+        variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(fastafile));
+        variables["[extension]"] = m->getExtension(fastafile);
+               string trimSeqFile = getOutputFileName("fasta",variables);
                outputNames.push_back(trimSeqFile); outputTypes["fasta"].push_back(trimSeqFile);
-        
-        string badSeqFile = thisOutputDir + m->getRootName(m->getSimpleName(fastafile)) + "pcr.scrap.fasta";
+        variables["[tag]"] = "scrap";
+        string badSeqFile = getOutputFileName("fasta",variables);
                
                
         length = 0;
-               if(oligosfile != ""){    readOligos();     }  if (m->control_pressed) {  return 0; }
+               if(oligosfile != ""){    readOligos();     if (m->debug) { m->mothurOut("[DEBUG]: read oligos file. numprimers = " + toString(primers.size()) + ", revprimers = " + toString(revPrimer.size()) + ".\n"); } }  if (m->control_pressed) {  return 0; }
         if(ecolifile != "") {    readEcoli();      }  if (m->control_pressed) {  return 0; }
         
         vector<unsigned long long> positions; 
@@ -313,7 +369,9 @@ int PcrSeqsCommand::execute(){
         if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {       m->mothurRemove(outputNames[i]); } return 0; }
                if (taxfile != "")                      {               readTax(badNames);              }
                if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        m->mothurRemove(outputNames[i]); } return 0; }
-        
+               if (countfile != "")                    {               readCount(badNames);            }
+               if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        m->mothurRemove(outputNames[i]); } return 0; }
+      
         m->mothurOutEndLine();
                m->mothurOut("Output File Names: "); m->mothurOutEndLine();
                for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
@@ -347,6 +405,11 @@ int PcrSeqsCommand::execute(){
                        if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTaxonomyFile(current); }
                }
         
+        itTypes = outputTypes.find("count");
+               if (itTypes != outputTypes.end()) {
+                       if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setCountTableFile(current); }
+               }
+        
                m->mothurOut("It took " + toString(time(NULL) - start) + " secs to screen " + toString(numFastaSeqs) + " sequences.");
                m->mothurOutEndLine();
 
@@ -442,7 +505,7 @@ int PcrSeqsCommand::createProcesses(string filename, string goodFileName, string
             if (i!=0) {extension += toString(i) + ".temp"; processIDS.push_back(i); }
             
                        // Allocate memory for thread data.
-                       pcrData* tempPcr = new pcrData(filename, goodFileName+extension, badFileName+extension, m, oligosfile, ecolifile, primers, revPrimer, nomatch, keepprimer, keepdots, start, end, length, lines[i].start, lines[i].end);
+                       pcrData* tempPcr = new pcrData(filename, goodFileName+extension, badFileName+extension, m, oligosfile, ecolifile, primers, revPrimer, nomatch, keepprimer, keepdots, start, end, length, pdiffs, lines[i].start, lines[i].end);
                        pDataArray.push_back(tempPcr);
                        
                        //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier
@@ -459,6 +522,9 @@ int PcrSeqsCommand::createProcesses(string filename, string goodFileName, string
                //Close all thread handles and free memory allocations.
                for(int i=0; i < pDataArray.size(); i++){
                        num += pDataArray[i]->count;
+            if (pDataArray[i]->count != pDataArray[i]->fend) {
+                m->mothurOut("[ERROR]: process " + toString(i) + " only processed " + toString(pDataArray[i]->count) + " of " + toString(pDataArray[i]->fend) + " sequences assigned to it, quitting. \n"); m->control_pressed = true; 
+            }
             for (set<string>::iterator it = pDataArray[i]->badSeqNames.begin(); it != pDataArray[i]->badSeqNames.end(); it++) {        badSeqNames.insert(*it);       }
                        CloseHandle(hThreadArray[i]);
                        delete pDataArray[i];
@@ -501,6 +567,10 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta
                int count = 0;
         set<int> lengths;
         
+        //pdiffs, bdiffs, primers, barcodes, revPrimers
+        map<string, int> faked;
+        TrimOligos trim(pdiffs, 0, primers, faked, revPrimer);
+        
                while (!done) {
             
                        if (m->control_pressed) {  break; }
@@ -510,6 +580,8 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta
             string trashCode = "";
                        if (currSeq.getName() != "") {
                 
+                if (m->debug) { m->mothurOut("[DEBUG]: seq name = " + currSeq.getName() + ".\n"); } 
+                
                 bool goodSeq = true;
                 if (oligosfile != "") {
                     map<int, int> mapAligned;
@@ -518,7 +590,7 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta
                     //process primers
                     if (primers.size() != 0) {
                         int primerStart = 0; int primerEnd = 0;
-                        bool good = findForward(currSeq, primerStart, primerEnd);
+                        bool good = trim.findForward(currSeq, primerStart, primerEnd);
                         
                         if(!good){     if (nomatch == "reject") { goodSeq = false; } trashCode += "f"; }
                         else{
@@ -530,8 +602,9 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta
                                 } 
                                 else                {  
                                     if (keepdots)   { currSeq.filterToPos(mapAligned[primerStart]);  }
-                                    else            { currSeq.setAligned(currSeq.getAligned().substr(mapAligned[primerStart]));                                              }
+                                    else            { currSeq.setAligned(currSeq.getAligned().substr(mapAligned[primerStart]));                                             }
                                 }
+                                isAligned(currSeq.getAligned(), mapAligned);
                             }else { 
                                 if (!keepprimer)    { currSeq.setAligned(currSeq.getUnaligned().substr(primerEnd)); } 
                                 else                { currSeq.setAligned(currSeq.getUnaligned().substr(primerStart)); } 
@@ -542,10 +615,10 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta
                     //process reverse primers
                     if (revPrimer.size() != 0) {
                         int primerStart = 0; int primerEnd = 0;
-                        bool good = findReverse(currSeq, primerStart, primerEnd);
+                        bool good = trim.findReverse(currSeq, primerStart, primerEnd);
                         if(!good){     if (nomatch == "reject") { goodSeq = false; } trashCode += "r"; }
                         else{ 
-                            //are you aligned
+                              //are you aligned
                             if (aligned) { 
                                 if (!keepprimer)    {  
                                     if (keepdots)   { currSeq.filterFromPos(mapAligned[primerStart]); }
@@ -640,75 +713,9 @@ int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta
        }
 }
 //********************************************************************/
-bool PcrSeqsCommand::findForward(Sequence& seq, int& primerStart, int& primerEnd){
-       try {
-               
-               string rawSequence = seq.getUnaligned();
-               
-               for(int j=0;j<primers.size();j++){
-                       string oligo = primers[j];
-                       
-                       if(rawSequence.length() < oligo.length()) {  break;  }
-                       
-                       //search for primer
-            int olength = oligo.length();
-            for (int j = 0; j < rawSequence.length()-olength; j++){
-                if (m->control_pressed) {  primerStart = 0; primerEnd = 0; return false; }
-                string rawChunk = rawSequence.substr(j, olength);
-                if(compareDNASeq(oligo, rawChunk)) {
-                    primerStart = j;
-                    primerEnd = primerStart + olength;
-                    return true;
-                }
-                
-            }
-        }      
-               
-        primerStart = 0; primerEnd = 0;
-               return false;
-               
-       }
-       catch(exception& e) {
-               m->errorOut(e, "TrimOligos", "stripForward");
-               exit(1);
-       }
-}
-//******************************************************************/
-bool PcrSeqsCommand::findReverse(Sequence& seq, int& primerStart, int& primerEnd){
-       try {
-               
-               string rawSequence = seq.getUnaligned();
-               
-               for(int i=0;i<revPrimer.size();i++){
-                       string oligo = revPrimer[i];
-                       if(rawSequence.length() < oligo.length()) {  break;  }
-                       
-                       //search for primer
-            int olength = oligo.length();
-            for (int j = rawSequence.length()-olength; j >= 0; j--){
-                 if (m->control_pressed) {  primerStart = 0; primerEnd = 0; return false; }
-                string rawChunk = rawSequence.substr(j, olength);
-            
-                if(compareDNASeq(oligo, rawChunk)) {
-                    primerStart = j;
-                    primerEnd = primerStart + olength;
-                    return true;
-                }
-                
-            }
-               }       
-               
-        primerStart = 0; primerEnd = 0;
-               return false;
-       }
-       catch(exception& e) {
-               m->errorOut(e, "PcrSeqsCommand", "findReverse");
-               exit(1);
-       }
-}
-//********************************************************************/
 bool PcrSeqsCommand::isAligned(string seq, map<int, int>& aligned){
        try {
+        aligned.clear();
         bool isAligned = false;
         
         int countBases = 0;
@@ -772,6 +779,7 @@ bool PcrSeqsCommand::readOligos(){
                m->openInputFile(oligosfile, inOligos);
                
                string type, oligo, group;
+        int primerCount = 0;
                
                while(!inOligos.eof()){
             
@@ -796,10 +804,10 @@ bool PcrSeqsCommand::readOligos(){
                                        // get rest of line in case there is a primer name
                                        while (!inOligos.eof()) {       
                         char c = inOligos.get(); 
-                        if (c == 10 || c == 13){       break;  } 
+                        if (c == 10 || c == 13 || c == -1){    break;  }
                         else if (c == 32 || c == 9){;} //space or tab
                                        } 
-                                       primers.push_back(oligo);
+                                       primers[oligo] = primerCount; primerCount++;
                 }else if(type == "REVERSE"){
                     string oligoRC = reverseOligo(oligo);
                     revPrimer.push_back(oligoRC);
@@ -854,7 +862,9 @@ int PcrSeqsCommand::writeAccnos(set<string> badNames){
        try {
                string thisOutputDir = outputDir;
                if (outputDir == "") {  thisOutputDir += m->hasPath(fastafile);  }
-               string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(fastafile)) + "bad.accnos";
+        map<string, string> variables; 
+               variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(fastafile));
+               string outputFileName = getOutputFileName("accnos",variables);
         outputNames.push_back(outputFileName); outputTypes["accnos"].push_back(outputFileName);
         
         ofstream out;
@@ -873,50 +883,16 @@ int PcrSeqsCommand::writeAccnos(set<string> badNames){
                exit(1);
        }
     
-}
-//******************************************************************/
-bool PcrSeqsCommand::compareDNASeq(string oligo, string seq){
-       try {
-               bool success = 1;
-               int length = oligo.length();
-               
-               for(int i=0;i<length;i++){
-                       
-                       if(oligo[i] != seq[i]){
-                               if(oligo[i] == 'A' || oligo[i] == 'T' || oligo[i] == 'G' || oligo[i] == 'C')    {       success = 0;    }
-                               else if((oligo[i] == 'N' || oligo[i] == 'I') && (seq[i] == 'N'))                                {       success = 0;    }
-                               else if(oligo[i] == 'R' && (seq[i] != 'A' && seq[i] != 'G'))                                    {       success = 0;    }
-                               else if(oligo[i] == 'Y' && (seq[i] != 'C' && seq[i] != 'T'))                                    {       success = 0;    }
-                               else if(oligo[i] == 'M' && (seq[i] != 'C' && seq[i] != 'A'))                                    {       success = 0;    }
-                               else if(oligo[i] == 'K' && (seq[i] != 'T' && seq[i] != 'G'))                                    {       success = 0;    }
-                               else if(oligo[i] == 'W' && (seq[i] != 'T' && seq[i] != 'A'))                                    {       success = 0;    }
-                               else if(oligo[i] == 'S' && (seq[i] != 'C' && seq[i] != 'G'))                                    {       success = 0;    }
-                               else if(oligo[i] == 'B' && (seq[i] != 'C' && seq[i] != 'T' && seq[i] != 'G'))   {       success = 0;    }
-                               else if(oligo[i] == 'D' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'G'))   {       success = 0;    }
-                               else if(oligo[i] == 'H' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'C'))   {       success = 0;    }
-                               else if(oligo[i] == 'V' && (seq[i] != 'A' && seq[i] != 'C' && seq[i] != 'G'))   {       success = 0;    }                       
-                               
-                               if(success == 0)        {       break;   }
-                       }
-                       else{
-                               success = 1;
-                       }
-               }
-               
-               return success;
-       }
-       catch(exception& e) {
-               m->errorOut(e, "PcrSeqsCommand", "compareDNASeq");
-               exit(1);
-       }
-       
 }
 //***************************************************************************************************************
 int PcrSeqsCommand::readName(set<string>& names){
        try {
                string thisOutputDir = outputDir;
                if (outputDir == "") {  thisOutputDir += m->hasPath(namefile);  }
-               string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(namefile)) + "pcr" + m->getExtension(namefile);
+               map<string, string> variables; 
+               variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(namefile));
+        variables["[extension]"] = m->getExtension(namefile);
+               string outputFileName = getOutputFileName("name", variables);
         
                ofstream out;
                m->openOutputFile(outputFileName, out);
@@ -974,7 +950,10 @@ int PcrSeqsCommand::readGroup(set<string> names){
        try {
                string thisOutputDir = outputDir;
                if (outputDir == "") {  thisOutputDir += m->hasPath(groupfile);  }
-               string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(groupfile)) + "pcr" + m->getExtension(groupfile);
+               map<string, string> variables; 
+               variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(groupfile));
+        variables["[extension]"] = m->getExtension(groupfile);
+               string outputFileName = getOutputFileName("group", variables);
                
                ofstream out;
                m->openOutputFile(outputFileName, out);
@@ -1021,7 +1000,11 @@ int PcrSeqsCommand::readTax(set<string> names){
        try {
                string thisOutputDir = outputDir;
                if (outputDir == "") {  thisOutputDir += m->hasPath(taxfile);  }
-               string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(taxfile)) + "pcr" + m->getExtension(taxfile);
+               map<string, string> variables; 
+               variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(taxfile));
+        variables["[extension]"] = m->getExtension(taxfile);
+               string outputFileName = getOutputFileName("taxonomy", variables);
+
                ofstream out;
                m->openOutputFile(outputFileName, out);
         
@@ -1061,6 +1044,67 @@ int PcrSeqsCommand::readTax(set<string> names){
                exit(1);
        }
 }
+//***************************************************************************************************************
+int PcrSeqsCommand::readCount(set<string> badSeqNames){
+       try {
+               ifstream in;
+               m->openInputFile(countfile, in);
+               set<string>::iterator it;
+               
+               map<string, string> variables; 
+               variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(countfile));
+        variables["[extension]"] = m->getExtension(countfile);
+               string goodCountFile = getOutputFileName("count", variables);
+
+        outputNames.push_back(goodCountFile);  outputTypes["count"].push_back(goodCountFile);
+               ofstream goodCountOut;  m->openOutputFile(goodCountFile, goodCountOut);
+               
+        string headers = m->getline(in); m->gobble(in);
+        goodCountOut << headers << endl;
+        
+        string name, rest; int thisTotal, removedCount; removedCount = 0;
+        bool wroteSomething = false;
+        while (!in.eof()) {
+            
+                       if (m->control_pressed) { goodCountOut.close(); in.close(); m->mothurRemove(goodCountFile); return 0; }
+            
+                       in >> name; m->gobble(in); 
+            in >> thisTotal; m->gobble(in);
+            rest = m->getline(in); m->gobble(in);
+            
+                       if (badSeqNames.count(name) != 0) { removedCount+=thisTotal; }
+                       else{
+                wroteSomething = true;
+                               goodCountOut << name << '\t' << thisTotal << '\t' << rest << endl;
+                       }
+               }
+               in.close();
+               goodCountOut.close();
+        
+        if (m->control_pressed) { m->mothurRemove(goodCountFile);   }
+        
+        if (wroteSomething == false) {  m->mothurOut("Your count file contains only sequences from the .accnos file."); m->mothurOutEndLine(); }
+        
+        //check for groups that have been eliminated
+        CountTable ct;
+        if (ct.testGroups(goodCountFile)) {
+            ct.readTable(goodCountFile);
+            ct.printTable(goodCountFile);
+        }
+               
+               if (m->control_pressed) { m->mothurRemove(goodCountFile);   }
+        
+        m->mothurOut("Removed " + toString(removedCount) + " sequences from your count file."); m->mothurOutEndLine();
+
+               
+               return 0;
+        
+       }
+       catch(exception& e) {
+               m->errorOut(e, "PcrSeqsCommand", "readCOunt");
+               exit(1);
+       }
+}
 /**************************************************************************************/