]> git.donarmstrong.com Git - mothur.git/blobdiff - trimseqscommand.cpp
added homova command
[mothur.git] / trimseqscommand.cpp
index fce4763dcb891ef56a57f24644b9574ffa51cd1a..7267b91da4d2ec193508579e8226c33d6758e9c7 100644 (file)
 #include "needlemanoverlap.hpp"
 
 //**********************************************************************************************************************
+
 vector<string> TrimSeqsCommand::getValidParameters(){  
        try {
                string Array[] =  {"fasta", "flip", "oligos", "maxambig", "maxhomop", "group","minlength", "maxlength", "qfile", 
-                                                                       "qthreshold", "qwindowaverage", "qstepsize", "qwindowsize", "qaverage", "rollaverage", "allfiles", "qtrim","tdiffs", "pdiffs", "bdiffs", "processors", "outputdir","inputdir"};
+                                                                       "qthreshold", "qwindowaverage", "qstepsize", "qwindowsize", "qaverage", "rollaverage",
+                                                                       "keepfirst", "removelast",
+                                                                       "allfiles", "qtrim","tdiffs", "pdiffs", "bdiffs", "processors", "outputdir","inputdir"};
                vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
                return myArray;
        }
@@ -23,7 +26,9 @@ vector<string> TrimSeqsCommand::getValidParameters(){
                exit(1);
        }
 }
+
 //**********************************************************************************************************************
+
 TrimSeqsCommand::TrimSeqsCommand(){    
        try {
                abort = true;
@@ -38,7 +43,9 @@ TrimSeqsCommand::TrimSeqsCommand(){
                exit(1);
        }
 }
+
 //**********************************************************************************************************************
+
 vector<string> TrimSeqsCommand::getRequiredParameters(){       
        try {
                string Array[] =  {"fasta"};
@@ -50,7 +57,9 @@ vector<string> TrimSeqsCommand::getRequiredParameters(){
                exit(1);
        }
 }
+
 //**********************************************************************************************************************
+
 vector<string> TrimSeqsCommand::getRequiredFiles(){    
        try {
                vector<string> myArray;
@@ -61,6 +70,7 @@ vector<string> TrimSeqsCommand::getRequiredFiles(){
                exit(1);
        }
 }
+
 //***************************************************************************************************************
 
 TrimSeqsCommand::TrimSeqsCommand(string option)  {
@@ -74,8 +84,10 @@ TrimSeqsCommand::TrimSeqsCommand(string option)  {
                
                else {
                        //valid paramters for this command
-                       string AlignArray[] =  {"fasta", "flip", "group","oligos", "maxambig", "maxhomop", "minlength", "maxlength", "qfile", 
-                                                                       "qthreshold", "qwindowaverage", "qstepsize", "qwindowsize", "qaverage", "rollaverage", "allfiles", "qtrim","tdiffs", "pdiffs", "bdiffs", "processors", "outputdir","inputdir"};
+                       string AlignArray[] =  {        "fasta", "flip", "oligos", "maxambig", "maxhomop", "group","minlength", "maxlength", "qfile", 
+                                                               "qthreshold", "qwindowaverage", "qstepsize", "qwindowsize", "qaverage", "rollaverage",
+                                                               "keepfirst", "removelast",
+                                                               "allfiles", "qtrim","tdiffs", "pdiffs", "bdiffs", "processors", "outputdir","inputdir"};
                        
                        vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
                        
@@ -195,7 +207,7 @@ TrimSeqsCommand::TrimSeqsCommand(string option)  {
                        temp = validParameter.validFile(parameters, "qthreshold", false);       if (temp == "not found") { temp = "0"; }
                        convert(temp, qThreshold);
                        
-                       temp = validParameter.validFile(parameters, "qtrim", false);            if (temp == "not found") { temp = "F"; }
+                       temp = validParameter.validFile(parameters, "qtrim", false);            if (temp == "not found") { temp = "t"; }
                        qtrim = m->isTrue(temp);
 
                        temp = validParameter.validFile(parameters, "rollaverage", false);      if (temp == "not found") { temp = "0"; }
@@ -204,14 +216,20 @@ TrimSeqsCommand::TrimSeqsCommand(string option)  {
                        temp = validParameter.validFile(parameters, "qwindowaverage", false);if (temp == "not found") { temp = "0"; }
                        convert(temp, qWindowAverage);
 
-                       temp = validParameter.validFile(parameters, "qwindowsize", false);      if (temp == "not found") { temp = "100"; }
+                       temp = validParameter.validFile(parameters, "qwindowsize", false);      if (temp == "not found") { temp = "50"; }
                        convert(temp, qWindowSize);
 
-                       temp = validParameter.validFile(parameters, "qstepsize", false);                if (temp == "not found") { temp = "10"; }
+                       temp = validParameter.validFile(parameters, "qstepsize", false);        if (temp == "not found") { temp = "1"; }
                        convert(temp, qWindowStep);
 
                        temp = validParameter.validFile(parameters, "qaverage", false);         if (temp == "not found") { temp = "0"; }
                        convert(temp, qAverage);
+
+                       temp = validParameter.validFile(parameters, "keepfirst", false);        if (temp == "not found") { temp = "0"; }
+                       convert(temp, keepFirst);
+
+                       temp = validParameter.validFile(parameters, "removelast", false);       if (temp == "not found") { temp = "0"; }
+                       convert(temp, removeLast);
                        
                        temp = validParameter.validFile(parameters, "allfiles", false);         if (temp == "not found") { temp = "F"; }
                        allFiles = m->isTrue(temp);
@@ -244,28 +262,36 @@ TrimSeqsCommand::TrimSeqsCommand(string option)  {
                exit(1);
        }
 }
+
 //**********************************************************************************************************************
 
 void TrimSeqsCommand::help(){
        try {
-               m->mothurOut("The trim.seqs command reads a fastaFile and creates .....\n");
-               m->mothurOut("The trim.seqs command parameters are fasta, flip, oligos, group, maxambig, maxhomop, minlength, maxlength, qfile, qthreshold, qaverage, diffs, qtrim and allfiles.\n");
+               m->mothurOut("The trim.seqs command reads a fastaFile and creates 2 new fasta files, .trim.fasta and scrap.fasta, as well as group files if you provide and oligos file.\n");
+               m->mothurOut("The .trim.fasta contains sequences that meet your requirements, and the .scrap.fasta contains those which don't.\n");
+               m->mothurOut("The trim.seqs command parameters are fasta, flip, oligos, group, maxambig, maxhomop, minlength, maxlength, qfile, qthreshold, qaverage, diffs, qtrim, keepfirst, removelast and allfiles.\n");
                m->mothurOut("The fasta parameter is required.\n");
                m->mothurOut("The group parameter allows you to enter a group file for your fasta file.\n");
                m->mothurOut("The flip parameter will output the reverse compliment of your trimmed sequence. The default is false.\n");
-               m->mothurOut("The oligos parameter .... The default is "".\n");
-               m->mothurOut("The maxambig parameter .... The default is -1.\n");
-               m->mothurOut("The maxhomop parameter .... The default is 0.\n");
-               m->mothurOut("The minlength parameter .... The default is 0.\n");
-               m->mothurOut("The maxlength parameter .... The default is 0.\n");
+               m->mothurOut("The oligos parameter allows you to provide an oligos file.\n");
+               m->mothurOut("The maxambig parameter allows you to set the maximum number of ambigious bases allowed. The default is -1.\n");
+               m->mothurOut("The maxhomop parameter allows you to set a maximum homopolymer length. \n");
+               m->mothurOut("The minlength parameter allows you to set and minimum sequence length. \n");
+               m->mothurOut("The maxlength parameter allows you to set and maximum sequence length. \n");
                m->mothurOut("The tdiffs parameter is used to specify the total number of differences allowed in the sequence. The default is pdiffs + bdiffs.\n");
                m->mothurOut("The bdiffs parameter is used to specify the number of differences allowed in the barcode. The default is 0.\n");
                m->mothurOut("The pdiffs parameter is used to specify the number of differences allowed in the primer. The default is 0.\n");
-               m->mothurOut("The qfile parameter .....\n");
-               m->mothurOut("The qthreshold parameter .... The default is 0.\n");
-               m->mothurOut("The qaverage parameter .... The default is 0.\n");
-               m->mothurOut("The allfiles parameter .... The default is F.\n");
-               m->mothurOut("The qtrim parameter .... The default is F.\n");
+               m->mothurOut("The qfile parameter allows you to provide a quality file.\n");
+               m->mothurOut("The qthreshold parameter allows you to set a minimum quality score allowed. \n");
+               m->mothurOut("The qaverage parameter allows you to set a minimum average quality score allowed. \n");
+               m->mothurOut("The qwindowsize parameter allows you to set a number of bases in a window. Default=50.\n");
+               m->mothurOut("The qwindowaverage parameter allows you to set a minimum average quality score allowed over a window. \n");
+               m->mothurOut("The rollaverage parameter allows you to set a minimum rolling average quality score allowed over a window. \n");
+               m->mothurOut("The qstepsize parameter allows you to set a number of bases to move the window over. Default=1.\n");
+               m->mothurOut("The allfiles parameter will create separate group and fasta file for each grouping. The default is F.\n");
+               m->mothurOut("The qtrim parameter will trim sequence from the point that they fall below the qthreshold and put it in the .trim file if set to true. The default is T.\n");
+               m->mothurOut("The keepfirst parameter trims the sequence to the first keepfirst number of bases after the barcode or primers are removed, before the sequence is checked to see if it meets the other requirements. \n");
+               m->mothurOut("The removelast removes the last removelast number of bases after the barcode or primers are removed, before the sequence is checked to see if it meets the other requirements.\n");
                m->mothurOut("The trim.seqs command should be in the following format: \n");
                m->mothurOut("trim.seqs(fasta=yourFastaFile, flip=yourFlip, oligos=yourOligos, maxambig=yourMaxambig,  \n");
                m->mothurOut("maxhomop=yourMaxhomop, minlength=youMinlength, maxlength=yourMaxlength)  \n");    
@@ -337,7 +363,7 @@ int TrimSeqsCommand::execute(){
                        outputNames.push_back(groupFile); outputTypes["group"].push_back(groupFile);
                        getOligos(fastaFileNames, qualFileNames);
                }
-               cout << fastaFileNames.size() << '\t' << qualFileNames.size() << endl;
+
                vector<unsigned long int> fastaFilePos;
                vector<unsigned long int> qFilePos;
                
@@ -360,18 +386,21 @@ int TrimSeqsCommand::execute(){
                #endif
                
                if (m->control_pressed) {  return 0; }                  
-               cout << "done with driver " << endl;                                                                    
-               for(int i=0;i<fastaFileNames.size();i++){
-                       cout << fastaFileNames[i] << endl;
                        
-                       if (m->isBlank(fastaFileNames[i])) { cout << fastaFileNames[i] << " was blank" << endl; remove(fastaFileNames[i].c_str());      }
-                       else if (filesToRemove.count(fastaFileNames[i]) > 0) { cout << fastaFileNames[i] << " was on the remove list" << endl; remove(fastaFileNames[i].c_str()); }
+               set<string> blanks;
+               for(int i=0;i<fastaFileNames.size();i++){
+                       if (m->isBlank(fastaFileNames[i])) {   blanks.insert(fastaFileNames[i]);        }
+                       else if (filesToRemove.count(fastaFileNames[i]) > 0) {  remove(fastaFileNames[i].c_str()); }
                        else {
                                ifstream inFASTA;
                                string seqName;
                                m->openInputFile(fastaFileNames[i], inFASTA);
                                ofstream outGroups;
                                string outGroupFilename = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[i])) + "groups";
+                               
+                               //if the fastafile is on the blanks list then the groups file should be as well
+                               if (blanks.count(fastaFileNames[i]) != 0) { blanks.insert(outGroupFilename); }
+                               
                                m->openOutputFile(outGroupFilename, outGroups);
                                outputNames.push_back(outGroupFilename); outputTypes["group"].push_back(outGroupFilename);  
                                
@@ -382,7 +411,7 @@ int TrimSeqsCommand::execute(){
                                                if(itCombo->second == i){       thisGroup = itCombo->first;     combos.erase(itCombo);  break;  }
                                        }
                                }else{ thisGroup = groupVector[i]; }
-                               cout << thisGroup << '\t' << i  << '\t' << comboStarts << endl; 
+                                       
                                while(!inFASTA.eof()){
                                        if(inFASTA.get() == '>'){
                                                inFASTA >> seqName;
@@ -394,32 +423,18 @@ int TrimSeqsCommand::execute(){
                                inFASTA.close();
                        }
                }
-       cout << "done with fastaFileNames " << endl;            
+               
+               for (set<string>::iterator itBlanks = blanks.begin(); itBlanks != blanks.end(); itBlanks++) {  remove((*(itBlanks)).c_str()); }
+               
+               blanks.clear();
                if(qFileName != ""){
                        for(int i=0;i<qualFileNames.size();i++){
-                               cout << qualFileNames[i] << endl;
-                               if (m->isBlank(qualFileNames[i])) { cout << qualFileNames[i] << " was blank" << endl; remove(qualFileNames[i].c_str()); }
-                               else if (filesToRemove.count(qualFileNames[i]) > 0) { cout << qualFileNames[i] << " was on the remove list" << endl; remove(qualFileNames[i].c_str()); }
-                               else {
-                                       ifstream inQual;
-                                       string seqName;
-                                       m->openInputFile(qualFileNames[i], inQual);
-//                                     ofstream outGroups;
-//                                     
-//                                     string thisGroup = "";
-//                                     if (i > comboStarts) {
-//                                             map<string, int>::iterator itCombo;
-//                                             for(itCombo=combos.begin();itCombo!=combos.end(); itCombo++){
-//                                                     if(itCombo->second == i){       thisGroup = itCombo->first;     combos.erase(itCombo);  break;  }
-//                                             }
-//                                     }
-//                                     else{ thisGroup = groupVector[i]; }
-                                       
-                                       inQual.close();
-                               }
+                               if (m->isBlank(qualFileNames[i])) {  blanks.insert(qualFileNames[i]);   }
+                               else if (filesToRemove.count(qualFileNames[i]) > 0) {  remove(qualFileNames[i].c_str()); }
                        }
                }
-               cout << "done with qualFileNames " << endl;
+               
+               for (set<string>::iterator itBlanks = blanks.begin(); itBlanks != blanks.end(); itBlanks++) {  remove((*(itBlanks)).c_str()); }
                
                if (m->control_pressed) { 
                        for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str()); }
@@ -447,7 +462,7 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string
        try {
                
                ofstream outFASTA;
-               int able = m->openOutputFile(trimFile, outFASTA);
+               m->openOutputFile(trimFile, outFASTA);
                
                ofstream scrapFASTA;
                m->openOutputFile(scrapFile, scrapFASTA);
@@ -460,8 +475,6 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string
                }
                
                ofstream outGroups;
-               //vector<ofstream*> fastaFileNames;
-               //vector<ofstream*> qualFileNames;
                
                if (oligoFile != "") {          
                        m->openOutputFile(groupFile, outGroups);   
@@ -474,6 +487,19 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string
                ifstream qFile;
                if(qFileName != "")     {       m->openInputFile(qFileName, qFile);     qFile.seekg(qline->start);  }
                
+               
+               for (int i = 0; i < fastaNames.size(); i++) { //clears old file
+                       ofstream temp;
+                       m->openOutputFile(fastaNames[i], temp);
+                       temp.close();
+               }
+               for (int i = 0; i < qualNames.size(); i++) { //clears old file
+                       ofstream temp;
+                       m->openOutputFile(qualNames[i], temp);
+                       temp.close();
+               }
+               
+                       
                bool done = false;
                int count = 0;
        
@@ -507,25 +533,12 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string
                                string trashCode = "";
                                int currentSeqsDiffs = 0;
 
-                               if(qFileName != ""){
-                                       if(qThreshold != 0)                     {       success = currQual.stripQualThreshold(currSeq, qThreshold);                     }
-                                       else if(qAverage != 0)          {       success = currQual.cullQualAverage(currSeq, qAverage);                          }
-                                       else if(qRollAverage != 0)      {       success = currQual.stripQualRollingAverage(currSeq, qRollAverage);      }
-                                       else if(qWindowAverage != 0){   success = currQual.stripQualWindowAverage(currSeq, qWindowStep, qWindowSize, qWindowAverage);   }
-
-                                       if (qtrim == 1 && (origSeq.length() != currSeq.getUnaligned().length())) { 
-                                               success = 0; //if you don't want to trim and the sequence does not meet quality requirements, move to scrap
-                                       }
-
-                                       if(!success)                            {       trashCode += 'q';       }
-                               }
-                       
                                if(barcodes.size() != 0){
                                        success = stripBarcode(currSeq, currQual, groupBar);
                                        if(success > bdiffs)            {       trashCode += 'b';       }
                                        else{ currentSeqsDiffs += success;  }
                                }
-
+                               
                                if(numFPrimers != 0){
                                        success = stripForward(currSeq, currQual, groupPrime);
                                        if(success > pdiffs)            {       trashCode += 'f';       }
@@ -533,11 +546,36 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string
                                }
                                
                                if (currentSeqsDiffs > tdiffs)  {       trashCode += 't';   }
-
+                               
                                if(numRPrimers != 0){
                                        success = stripReverse(currSeq, currQual);
                                        if(!success)                            {       trashCode += 'r';       }
                                }
+
+                               if(keepFirst != 0){
+                                       success = keepFirstTrim(currSeq, currQual);
+                               }
+                               
+                               if(removeLast != 0){
+                                       success = removeLastTrim(currSeq, currQual);
+                                       if(!success)                            {       trashCode += 'l';       }
+                               }
+
+                               
+                               if(qFileName != ""){
+                                       int origLength = currSeq.getNumBases();
+                                       
+                                       if(qThreshold != 0)                     {       success = currQual.stripQualThreshold(currSeq, qThreshold);                     }
+                                       else if(qAverage != 0)          {       success = currQual.cullQualAverage(currSeq, qAverage);                          }
+                                       else if(qRollAverage != 0)      {       success = currQual.stripQualRollingAverage(currSeq, qRollAverage);      }
+                                       else if(qWindowAverage != 0){   success = currQual.stripQualWindowAverage(currSeq, qWindowStep, qWindowSize, qWindowAverage);   }
+                                       else                                            {       success = 1;                            }
+                                       
+                                       //you don't want to trim, if it fails above then scrap it
+                                       if ((!qtrim) && (origLength != currSeq.getNumBases())) {  success = 0; }
+                                       
+                                       if(!success)                            {       trashCode += 'q';       }
+                               }                               
                
                                if(minLength > 0 || maxLength > 0){
                                        success = cullLength(currSeq);
@@ -552,7 +590,10 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string
                                        if(!success)                            {       trashCode += 'n';       }
                                }
                                
-                               if(flip){       currSeq.reverseComplement();            currQual.flipQScores(); }               // should go last                       
+                               if(flip){               // should go last                       
+                                       currSeq.reverseComplement();
+                                       currQual.flipQScores(); 
+                               }
                                
                                if(trashCode.length() == 0){
                                        currSeq.setAligned(currSeq.getUnaligned());
@@ -718,7 +759,8 @@ int TrimSeqsCommand::createProcessesCreateTrim(string filename, string qFileName
                
                //append files
                for(int i=0;i<processIDS.size();i++){
-                       m->mothurOut("Appending files from process " + processIDS[i]); m->mothurOutEndLine();
+                       
+                       m->mothurOut("Appending files from process " + toString(processIDS[i])); m->mothurOutEndLine();
                        
                        m->appendFiles((trimFile + toString(processIDS[i]) + ".temp"), trimFile);
                        remove((trimFile + toString(processIDS[i]) + ".temp").c_str());
@@ -753,7 +795,7 @@ int TrimSeqsCommand::createProcessesCreateTrim(string filename, string qFileName
                                }
                        }
                        
-                       if (allFiles) { m->mothurOut("Done with allfile"); m->mothurOutEndLine(); }
+                       if (allFiles) { m->mothurOut("Done with allfiles"); m->mothurOutEndLine(); }
                }
        
                return exitCommand;
@@ -848,6 +890,7 @@ int TrimSeqsCommand::setLines(string filename, string qfilename, vector<unsigned
                exit(1);
        }
 }
+
 //***************************************************************************************************************
 
 void TrimSeqsCommand::getOligos(vector<string>& outFASTAVec, vector<string>& outQualVec){
@@ -862,7 +905,7 @@ void TrimSeqsCommand::getOligos(vector<string>& outFASTAVec, vector<string>& out
                //int indexPrimer = 0;
                
                while(!inOligos.eof()){
-                       inOligos >> type;
+                       inOligos >> type; m->gobble(inOligos);
                                        
                        if(type[0] == '#'){
                                while (!inOligos.eof()) {       char c = inOligos.get(); if (c == 10 || c == 13){       break;  }       } // get rest of line if there's any crap there
@@ -893,29 +936,29 @@ void TrimSeqsCommand::getOligos(vector<string>& outFASTAVec, vector<string>& out
                                        map<string, int>::iterator itPrime = primers.find(oligo);
                                        if (itPrime != primers.end()) { m->mothurOut("primer " + oligo + " is in your oligos file already."); m->mothurOutEndLine();  }
                                        
-                                       primers[oligo]=index; index++;
-                                       groupVector.push_back(group);
+                                               primers[oligo]=index; index++;
+                                               groupVector.push_back(group);
                                        
-                                       if(allFiles){
-                                               outFASTAVec.push_back((outputDir + m->getRootName(m->getSimpleName(fastaFile)) + group + ".fasta"));
-                                               if(qFileName != ""){
-                                                       outQualVec.push_back((outputDir + m->getRootName(m->getSimpleName(qFileName)) + group + ".qual"));
-                                               }
-                                               if (group == "") { //if there is not a group for this primer, then this file will not get written to, but we add it to keep the indexes correct
-                                                       filesToRemove.insert((outputDir + m->getRootName(m->getSimpleName(fastaFile)) + group + ".fasta"));
+                                               if(allFiles){
+                                                       outFASTAVec.push_back((outputDir + m->getRootName(m->getSimpleName(fastaFile)) + group + ".fasta"));
                                                        if(qFileName != ""){
-                                                               filesToRemove.insert((outputDir + m->getRootName(m->getSimpleName(qFileName)) + group + ".qual"));
+                                                               outQualVec.push_back((outputDir + m->getRootName(m->getSimpleName(qFileName)) + group + ".qual"));
+                                                       }
+                                                       if (group == "") { //if there is not a group for this primer, then this file will not get written to, but we add it to keep the indexes correct
+                                                               filesToRemove.insert((outputDir + m->getRootName(m->getSimpleName(fastaFile)) + group + ".fasta"));
+                                                               if(qFileName != ""){
+                                                                       filesToRemove.insert((outputDir + m->getRootName(m->getSimpleName(qFileName)) + group + ".qual"));
+                                                               }
+                                                       }else {
+                                                               outputNames.push_back((outputDir + m->getRootName(m->getSimpleName(fastaFile)) + group + ".fasta"));
+                                                               outputTypes["fasta"].push_back((outputDir + m->getRootName(m->getSimpleName(fastaFile)) + group + ".fasta"));
+                                                               if(qFileName != ""){
+                                                                       outputNames.push_back((outputDir + m->getRootName(m->getSimpleName(qFileName)) + group + ".qual"));
+                                                                       outputTypes["qual"].push_back((outputDir + m->getRootName(m->getSimpleName(qFileName)) + group + ".qual"));
+                                                               }                                                       
                                                        }
-                                               }else {
-                                                       outputNames.push_back((outputDir + m->getRootName(m->getSimpleName(fastaFile)) + group + ".fasta"));
-                                                       outputTypes["fasta"].push_back((outputDir + m->getRootName(m->getSimpleName(fastaFile)) + group + ".fasta"));
-                                                       if(qFileName != ""){
-                                                               outputNames.push_back((outputDir + m->getRootName(m->getSimpleName(qFileName)) + group + ".qual"));
-                                                               outputTypes["qual"].push_back((outputDir + m->getRootName(m->getSimpleName(qFileName)) + group + ".qual"));
-                                                       }                                                       
                                                }
-                                       }
-
+                                       
                                }
                                else if(type == "REVERSE"){
                                        Sequence oligoRC("reverse", oligo);
@@ -929,19 +972,20 @@ void TrimSeqsCommand::getOligos(vector<string>& outFASTAVec, vector<string>& out
                                        map<string, int>::iterator itBar = barcodes.find(oligo);
                                        if (itBar != barcodes.end()) { m->mothurOut("barcode " + oligo + " is in your oligos file already."); m->mothurOutEndLine();  }
                                        
-                                       barcodes[oligo]=index; index++;
-                                       groupVector.push_back(group);
+                                               barcodes[oligo]=index; index++;
+                                               groupVector.push_back(group);
+                                               
+                                               if(allFiles){
+                                                       outputNames.push_back((outputDir + m->getRootName(m->getSimpleName(fastaFile)) + group + ".fasta"));
+                                                       outputNames.push_back((outputDir + m->getRootName(m->getSimpleName(fastaFile)) + group + ".fasta"));
+                                                       outFASTAVec.push_back((outputDir + m->getRootName(m->getSimpleName(fastaFile)) + group + ".fasta"));
+                                                       if(qFileName != ""){
+                                                               outQualVec.push_back((outputDir + m->getRootName(m->getSimpleName(qFileName)) + group + ".qual"));
+                                                               outputNames.push_back((outputDir + m->getRootName(m->getSimpleName(qFileName)) + group + ".qual"));
+                                                               outputTypes["qual"].push_back((outputDir + m->getRootName(m->getSimpleName(qFileName)) + group + ".qual"));
+                                                       }                                                       
+                                               }
                                        
-                                       if(allFiles){
-                                               outputNames.push_back((outputDir + m->getRootName(m->getSimpleName(fastaFile)) + group + ".fasta"));
-                                               outputNames.push_back((outputDir + m->getRootName(m->getSimpleName(fastaFile)) + group + ".fasta"));
-                                               outFASTAVec.push_back((outputDir + m->getRootName(m->getSimpleName(fastaFile)) + group + ".fasta"));
-                                               if(qFileName != ""){
-                                                       outQualVec.push_back((outputDir + m->getRootName(m->getSimpleName(qFileName)) + group + ".qual"));
-                                                       outputNames.push_back((outputDir + m->getRootName(m->getSimpleName(qFileName)) + group + ".qual"));
-                                                       outputTypes["qual"].push_back((outputDir + m->getRootName(m->getSimpleName(qFileName)) + group + ".qual"));
-                                               }                                                       
-                                       }
                                }else{  m->mothurOut(type + " is not recognized as a valid type. Choices are forward, reverse, and barcode. Ignoring " + oligo + "."); m->mothurOutEndLine();  }
                        }
                        m->gobble(inOligos);
@@ -1058,7 +1102,6 @@ int TrimSeqsCommand::stripBarcode(Sequence& seq, QualityScores& qual, int& group
                                oligo = oligo.substr(0,alnLength);
                                temp = temp.substr(0,alnLength);
                                
-                               int newStart=0;
                                int numDiff = countDiffs(oligo, temp);
                                
 //                             cout << oligo << '\t' << temp << '\t' << numDiff << endl;                               
@@ -1127,7 +1170,6 @@ int TrimSeqsCommand::stripForward(Sequence& seq, QualityScores& qual, int& group
                                seq.setUnaligned(rawSequence.substr(oligo.length()));
                                if(qual.getName() != ""){
                                        qual.trimQScores(oligo.length(), -1);
-                                       
                                }
                                success = 0;
                                break;
@@ -1184,7 +1226,6 @@ int TrimSeqsCommand::stripForward(Sequence& seq, QualityScores& qual, int& group
                                oligo = oligo.substr(0,alnLength);
                                temp = temp.substr(0,alnLength);
                                
-                               int newStart=0;
                                int numDiff = countDiffs(oligo, temp);
                                
 //                             cout << oligo << '\t' << temp << '\t' << numDiff << endl;                               
@@ -1265,6 +1306,52 @@ bool TrimSeqsCommand::stripReverse(Sequence& seq, QualityScores& qual){
 
 //***************************************************************************************************************
 
+bool TrimSeqsCommand::keepFirstTrim(Sequence& sequence, QualityScores& qscores){
+       try {
+               bool success = 1;
+               if(qscores.getName() != ""){
+                       qscores.trimQScores(-1, keepFirst);
+               }
+               sequence.trim(keepFirst);
+               return success;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "keepFirstTrim", "countDiffs");
+               exit(1);
+       }
+       
+}      
+
+//***************************************************************************************************************
+
+bool TrimSeqsCommand::removeLastTrim(Sequence& sequence, QualityScores& qscores){
+       try {
+               bool success = 0;
+               
+               int length = sequence.getNumBases() - removeLast;
+               
+               if(length > 0){
+                       if(qscores.getName() != ""){
+                               qscores.trimQScores(-1, length);
+                       }
+                       sequence.trim(length);
+                       success = 1;
+               }
+               else{
+                       success = 0;
+               }
+
+               return success;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "removeLastTrim", "countDiffs");
+               exit(1);
+       }
+       
+}      
+
+//***************************************************************************************************************
+
 bool TrimSeqsCommand::cullLength(Sequence& seq){
        try {
        
@@ -1361,6 +1448,7 @@ bool TrimSeqsCommand::compareDNASeq(string oligo, string seq){
        }
 
 }
+
 //***************************************************************************************************************
 
 int TrimSeqsCommand::countDiffs(string oligo, string seq){
@@ -1396,78 +1484,5 @@ int TrimSeqsCommand::countDiffs(string oligo, string seq){
        }
 
 }
-//***************************************************************************************************************
-
-//bool TrimSeqsCommand::stripQualThreshold(Sequence& seq, ifstream& qFile){
-//     try {
-//             
-//             string rawSequence = seq.getUnaligned();
-//             int seqLength = seq.getNumBases();
-//             bool success = 0;       //guilty until proven innocent
-//             string name;
-//             
-//             qFile >> name;
-//             if (name[0] == '>') {  if(name.substr(1) != seq.getName())      {       m->mothurOut("sequence name mismatch btwn fasta: " + seq.getName() + " and qual file: " + name); m->mothurOutEndLine(); }  }
-//             
-//             while (!qFile.eof())    {       char c = qFile.get(); if (c == 10 || c == 13){  break;  }       }
-//             
-//             int score;
-//             int end = seqLength;
-//             
-//             for(int i=0;i<seqLength;i++){
-//                     qFile >> score;
-//                     
-//                     if(score < qThreshold){
-//                             end = i;
-//                             break;
-//                     }
-//             }
-//             for(int i=end+1;i<seqLength;i++){
-//                     qFile >> score;
-//             }
-//             
-//             seq.setUnaligned(rawSequence.substr(0,end));
-//             
-//             return 1;
-//     }
-//     catch(exception& e) {
-//             m->errorOut(e, "TrimSeqsCommand", "stripQualThreshold");
-//             exit(1);
-//     }
-//}
-
-//***************************************************************************************************************
-
-//bool TrimSeqsCommand::cullQualAverage(Sequence& seq, ifstream& qFile){
-//     try {
-//             string rawSequence = seq.getUnaligned();
-//             int seqLength = seq.getNumBases();
-//             bool success = 0;       //guilty until proven innocent
-//             string name;
-//             
-//             qFile >> name;
-//             if (name[0] == '>') {  if(name.substr(1) != seq.getName())      {       m->mothurOut("sequence name mismatch btwn fasta: " + seq.getName() + " and qual file: " + name); m->mothurOutEndLine(); }  }
-//             
-//             while (!qFile.eof())    {       char c = qFile.get(); if (c == 10 || c == 13){  break;  }       }
-//             
-//             float score;    
-//             float average = 0;
-//             
-//             for(int i=0;i<seqLength;i++){
-//                     qFile >> score;
-//                     average += score;
-//             }
-//             average /= seqLength;
-//
-//             if(average >= qAverage) {       success = 1;    }
-//             else                                    {       success = 0;    }
-//             
-//             return success;
-//     }
-//     catch(exception& e) {
-//             m->errorOut(e, "TrimSeqsCommand", "cullQualAverage");
-//             exit(1);
-//     }
-//}
 
 //***************************************************************************************************************